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ABSTRACT 


Transformation optics is the current method used to design metamaterial structures that 
manipulate the path of electromagnetic radiation. This approach, however, relies upon 
a completely linear response of the polarization and magnetization fields with respect to 
incident electromagnetic field intensities. As those field intensities rise, such as from a 
hypothetical directed energy weapon, nonlinear effects, which are unaccounted for in a 
completely linear theory, are observed. In order to investigate the behavior of a transforma¬ 
tion optics-derived structure in such a high-field intensity regime, we propose to employ 
an iterative solution to the Maxwell equations for such a structure, and compare these 
results to those of the purely linear transformation optics model. Examining the first-order 
results of this approach, we observe a strong dependence of response field amplitude 
upon the wavelength of incident radiation. 
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CHAPTER 1: 
Introduction 


As the Department of Defense (DOD) continues its development and implementation of 
a new category of weaponry, directed energy weapons (DEW), other nations are expected 
to follow suit [1]. The design method known as transformation optics (TO), developed by 
Pendry, Schurig, and Smith [2], provides a highly useful design tool for determining the 
electromagnetic material properties required in order to redirect incoming electromagnetic 
radiation along a more desirable path. This design method has proven effective in the 
fabrication of real world structures using metamaterials — materials engineered to 
possess properties which are not found in nature [3]. Although both powerful and 
elegant in its formulation, TO relies upon a purely linear formulation of the response 
of materials to electric and magnetic fields. It has been known since at least 1941 [4] 
that sufficiently large electric and magnetic field amplitudes induce nonlinear (i.e., second 
and higher order) material responses. As these high field amplitudes are naturally expected 
in an environment in which the concern is the countering of certain classes of DEW, 
such as high power lasers and microwave devices, a question of the effectiveness of 
TO derived redirection structures for counter-directed energy weapon (CDEW) 
applications arises. We attempt here to address this question of the utility of the 
transformation optics approach in a CDEW environment. 


1.1 Directed Energy Weapons 

A directed energy weapon is a weapon which delivers electromagnetic energy, rather than 
a traditional projectile, to its target. Typical examples of directed energy weapons include 
high power microwave (HPM) or high power radio frequency (HERE) devices, designed to 
induce large voltages and currents in adversary electronic systems in order to render those 
systems inoperative, and high energy lasers (HEE), which ablate target material through 
the delivery of a focused high intensity beam of electromagnetic radiation onto a small 
area of the target. These weapon systems have been active areas of Navy research 
since the 1960s, and are the subject of ongoing development [1], [5]. As research into 
these devices proceeds, concern has grown regarding the ability of potential adversaries 
to field similar 
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weaponry, particularly within the realm of anti satellite applications [6]. 


1.2 Counter-Directed Energy Weapons 

In response to this concern, the Office of Naval Research (ONR), in conjunction with the 
Naval Postgraduate School (NPS), United States Naval Academy (USNA), and Naval Re¬ 
search Laboratory (NRL) is "investigating basic research topics in counter threats from 
directed energy weapons systems, such as high-power lasers or microwaves" [7]. These 
research topics include, from [1], 

a. Advanced materials including nano- and/or nonlinear materials for enhanced 
HEL protection of sensors, optics, airframe, etc. 

b. Metamaterial structures for the control and mitigation of HEL and HPRE 
irradiation. 

c. Techniques for HEL mitigation such as use of plasmas and obscurants. 

d. HEL protection by degrading atmospheric transmission (e.g. thermal 
blooming, scattering, absorption aids, and turbulence). 

e. Modeling and sensing of off-axis detection and source geo-location. 

f. Novel instrumentation for detection of HEL and HPRE irradiation. 

g. Active/passive circuit protection and limiters for HPRE 

h. Modeling of HPRE and HEL effects to materials, electronics and sensors as 
applied to CDEW objectives. 


1.2.1 Transformation Optics 

The field of transformation optics, founded in [2] by Pendry, Schurig, and Smith, uses the 
form invariance of Maxwell’s equations, the fundamental equations of electromagnetism, 
as well as the transformation properties of various electromagnetic fields and material prop¬ 
erties under coordinate transformations, in order to determine the material properties nec¬ 
essary to force electromagnetic fields to behave in some desired fashion. A design tool for 
directing electromagnetic radiation along a specified path, transformation optics’ potential 
for CDEW applications was quickly recognized. Just as quickly, the material parameters 
necessary to realize a working electromagnetic redirection structure (as seen below), were 
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discovered to be unlike those found in nature. Speeifieally, negative eleetrie and magnetie 
susceptibility values are commonly required [2], [3], [8], [9], [10]. In addition, several 
applications require unrealistieally diverging permittivity and permeability as an outer or 
inner edge of the strueture is approaehed, as well as eontinuously (spatially) varying permit¬ 
tivity and permeability [8], [9]. Figure 1.1 illustrates both ideal and approximate material 
parameters required to aehieve the cylindrical redirection structure of [9], as well as the 
simulated performanee of the same strueture with ideal material parameters. 



Figure 1.1: Left: Required ideal material properties for a cylindrical cloak of inner radius Ri 
and outer radius of 7?2 = 2Ri, as well as design compromises in order to make the structure 
realizeable. Right: Full-wave simulation of the performance of the idealized material parameters. 
Note that the waveform is undisturbed upon exiting the structure. Grey lines indicate the local 
direction of power flow. Source: [11]. 


While approximations and informed design eompromises may effeetively address the issue 
of diverging or continuously varying material properties [3], [12], the problem of material 
properties which are not observed in nature may be resolved through the use of metamate¬ 
rials [12], [13]. Of partieular note, however is the eompletely linear eonstitutive relation 


P‘=£oxi^Ej ( 1 . 1 ) 

= ( 1 . 2 ) 

employed in TO, where P is the polarization, x'e are the eomponents of the eleetrie sus- 
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ceptibility, E is the electric field, M is the magnetization, are the components of the 
magnetic susceptibility, and H is the magnetic field, and the Einstein summation conven¬ 
tion is in use. A treatment of the theory of transformation optics is provided in Chapter 3 of 
this work. 

1.2.2 Metamaterials 

A metamaterial is a structure which is designed to possess aggregate material 
properties which have not yet been observed in nature. Within the context of 
electromagnetic interac-tions, these novel material responses are typically designed 
through the addition of struc-tural features which are smaller than the wavelength of 
radiation to be scattered. Specific to the negative permittivity and permeability often 
required for implementation of TO applications, an array of split-ring resonators, structures 
featuring two conducting, concentric, incomplete rings, is commonly used. In 2009, a 
cylindrical cloak operating at microwave frequencies, designed through TO and employing 
an array of split-ring resonators among other design features, was successfully 
demonstrated [3]. See Figure 1.2 for an example of this design. 



4.7 4.8 4.9 5 5.1 5.2 5.3 5.4 5.5 


Frequency (GHz) 

Figure 1.2: Left: Diagram of a split-ring resonator alongside its resonance curve (source: [14]). 
Right: Cylindrical microwave frequency cloak (source: [3]), with diagrams of split-ring resonator 
design based on radial position, and diagram of idealized material parameters. 
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Metamaterials may therefore be viewed as a promising set of struetures with whieh we may 
realize a CDEW foeused eleetromagnetie redireetion deviee that has been designed through 
transformation opties. 


1.3 Nonlinear Optics 

As mentioned in the previous seetion, the TO design approaeh assumes a linear set of 
eonstitutive relations as given by (1.1) and (1.2). High amplitude eleetrie and magnetie 
fields, however, are known to induee seeond order and higher effeets typieally modeled by, 
as an extension of (1.1), 


= Co [tjEj + tJ^EjEk + xy^'EjEkEi +...] (1.3) 

where Xe^ the eomponents of the seeond-order eleetrie suseeptibility and Xe^^ are the 
eomponents of the third-order eleetrie suseeptibility [15]. Seeond-order effeets inelude fre- 
queney doubling, while third order effeets inelude the Kerr effeet, wherein refraetive index 
varies with field strength [15]. The linear eonstitutive relations (1.1) and (1.2) do not allow 
for explieit realization of these observed phenomena. Therefore, the TO approaeh does not 
take nonlinear optieal effeets into aeeount, and thus struetures designed through TO may 
be expeeted to stray from design expeetations as field strengths rise high enough for non¬ 
linear effeets to beeome non-negligible. As a potential DEW attaek upon a TO designed 
redireetion strueture may involve eleetrie and magnetie fields of suffieient intensity to trig¬ 
ger signifieant nonlinear effeets, the question of applieability of the TO design method for 
CDEW applieations arises. This question is the foeus of the following ehapters. 
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CHAPTER 2: 

Classical Electrodynamics in Matter 


As TO describes a design method through which we may control the behavior of elec¬ 
tromagnetic fields in matter, a description of the physical phenomenon of 
electromagnetic waves is in order. The objective of the following chapter will be the 
description of gen-eral electromagnetic phenomena via the Maxwell equations, followed 
by the application of these concepts to electromagnetic waves in free space as well as in 
media. 


2.1 Maxwell’s Equations 

Maxwell’s equations, the governing differential equations for all classical electrodynamic 
phenomena, are expressed as: 


WE = — 

£o 

(2.1) 

V. - 

V xE = —=- 
dt 

(2.2) 

WB = d 

(2.3) 

dE 

W X B = HoJ + £oHo— 

(2.4) 


where E is the electric field, B the magnetic flux density, p the charge density, J the cur¬ 
rent density, £„ the permittivity of vacuum, and po the permeability of vacuum. Note that 
overall charge conservation is consistent with the Maxwell equations as written: Taking the 
divergence of the Ampere-Maxwell law (2.4), 


V- (^Vx5j =poV-7 + eoAtoV--^ 


Since the divergence of a curl is always equal to zero, and the spatial and time derivatives 
are independent of each other. 
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diV-E 

0 = fiQV ■ J + £oAto 


( 2 . 6 ) 


Inserting Gauss’ law (2.1) into the above relation, 


£0 



dt 


-VJ 


dp 

dt 


-VJ 


(2.7) 


( 2 . 8 ) 


2.2 Maxwell’s Equations in Matter 

While the Maxwell equations (2.1)- (2.4) provide a valid deseription of all eleetromagnetie 
phenomena, they are somewhat eumbersome when eonsidering these fields in matter, where 
eharge may be bound up in a material’s eonstituent atoms. A more eonvenient approach 
is to first draw a distinction between py, the free charge density, and pt,, the bound charge 
density. Defining the polarization as 


(2.9) 

where pi is the electric dipole moment of a particular atom or molecule of interest. The 
bound charge density is defined as 


Pb = -VP 


( 2 . 10 ) 


If charge is to be conserved, then 


dpb 

dt 


-V-fp 


( 2 . 11 ) 
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where Jp is the polarization eurrent assoeiated with the motion of bound eharges. The 
above relationship thus requires 


Jn 


dP 

~di 


yielding 


dpb 

dt 




= -y-jp 


( 2 . 12 ) 


(2.13) 


The free eharge density pf is defined as any eharge density whieh is not part of the bound 
eharge density, giving 


pf + pb = p (2.14) 

Similarly, the magnetization M is defined 

M=^'£mi (2.15) 

where Mi is the magnetie moment of an atom or moleeule of interest. The bound eurrent 
density Jb is 

Jb = VxM (2.16) 

The free eurrent density Jf is the eurrent density unaeeounted for by the bound or polariza¬ 
tion eurrent densities. That is. 


J — J f Jp -\- Jb 


(2.17) 
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Inserting (2.14) into (2.1), 


y.E = 

£o 


(2.18) 


VE 


Pf 

£o £o 


(2.19) 


V.(eo£ + p) =pf (2.20) 

The quantity EqE + .P is the eleetrie displaeement field, denoted D. Arriving at Gauss’ Law 
for maeroseopie media. 


W b = pf 


Inserting (2.17) into (2.4), 


^ , , dE 

V X 5 — /io [7/ + 7p + 7fc j +eo/io-^ 


dP 


dE 


VxB = po\Jf + — + VxM]+ £oPo-^ 


Vx I- M =Jf + 

Po 


dt 


, B do 

Vx (-M =Jf + — 

Po 


dt 


( 2 . 21 ) 


( 2 . 22 ) 


(2.23) 


(2.24) 


(2.25) 
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B ^ 

The quantity-A/ is called the magnetic field and denoted H. The Ampere-Maxwell 

Bo 

law for macroscopic media is thus 


VxH = Jf + — (2.26) 

Replacing (2.1) with (2.21), and (2.4) with (2.26), we arrive at Maxwell’s equations in 
matter 


V b = pf 


(2.27) 




dB 


(2.28) 


Vfi = 0 


(2.29) 


VxH 



(2.30) 


2.3 The Wave Equations 

2.3.1 The Wave Equations in Vacuum 

Let us consider the Maxwell equations in vacuum in the absence of sources. Setting p = 0 
and 7 = 0, 


y-E = o 


(2.31) 
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dB 

~dt 


(2.32) 


VxE = - 


V-5 = 0 


(2.33) 


Vx5 


dE 


(2.34) 


Taking the curl of (2.32), 


V X V X £ 



(2.35) 


Employing the vector identity V x V x V 


_v2y + v(^V-y 


V^E + V V-E = 


(9 ( V X 5 


dt 


(2.36) 


Inserting (2.31) into the above relationship, 




d(vxB^ 

dt 


(2.37) 


Inserting the sourceless Ampere-Maxwell law (2.34) in vacuum into the above relationship. 


= 



(2.38) 


d^E 

~W 


Co Mo 


(2.39) 
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Taking the curl of (2.34), 


V X V X 5 


V 


dE 

xeoMo^ 


(2.40) 


-V^B + V 



d(vxE^ 


(2.41) 


Inserting Faraday’s law (2.32) into the above relationship, and noting from (2.33) that V ■ 
5 = 0, 





(2.42) 


d^B 


-^V^B 
Co Mo 


(2.43) 


Recognizing (2.39) and (2.43) as wave equations for electric and magnetic fields which 

each propagate with velocity ——= = c, which is experimentally equal to the speed of 

a/SoMo 

light in vacuum, we arrive at 


d^E 




(2.44) 


d^B 

w 




(2.45) 


2.3.2 The Wave Equations in Matter 

Taking E and H as our primary fields. Gauss’ law in matter (2.27) becomes 
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V • ( £qE +p] = Pf 


(2.46) 


^_?-V vP 

£o £o 


(2.47) 


From Faraday’s law in matter (2.28), 


V X P = -po- 


diH + M 


dt 


(2.48) 


Equation (2.29) gives 


AtoV- [H + M] =0 


(2.49) 


sj.H = -VM 


(2.50) 


The Ampere-Maxwell law in matter (2.30) may be expressed 


VxH=Jf + 


d ( £{)E + P 


dt 


(2.51) 


As was the case in vacuum, we assume our region of interest is free of sources - in this case, 
the free charge density and free current density pf and Jf, respectively. Taking pf = 0 and 
Jf = 0, Maxwell’s equations now reduce to 


VP = -—VP 
£0 


(2.52) 
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(2.53) 


V X £■ = -jiQ- 


dlH + M 


VH = -VM 


(2.54) 


Vx// = 


d ( EqE + P 


(2.55) 


Once again similar to the case in vaeuum, we take the eurl of the Maxwell equations in¬ 
volving the eurl of the primary fields. Taking the curl of (2.48), 


V xV X E — — /ioV X 


dlH+M 


(2.56) 


Again employing the vector identity V X V X y = — V^y -|- V ■ y j, and exploiting the 

independence of the temporal and spatial partial derivatives, 


- V^E -h V (V ■ £) = -Ho 


d Vx H + M 


(2.57) 


Inserting (2.52) and (2.55) into the above relationship. 


2 - 1 / -3 dleo£^+Pj 

-V^E - v(v-P)~-jUo^ —^- - + VxM 

Co V / dt 


(2.58) 


^2- 1^/^ d^E d^P / dM 


(2.59) 
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1 


Forming the wave operator on the left, and identifying c = 




as previously, 


^ d^E d^P dM\ 1 ^/^ - 

V^E-—— =^o—+ ^o V(V-F 


dt^ dt^ 


dt Co V 


(2.60) 


Taking the eurl of (2.55), 


VxVx// = Vx 


d ( £qE + P 


dt 


(2.61) 


-y^H+y [V-H] =£o 


d (V xE) ^fVxF 

+ ■ 


dt 


dt 


(2.62) 


Inserting (2.53) and (2.54) into the above relationship. 


V^H-V [V -M) = Co¬ 
ot 


-Po- 


dlH + M 


dt 


+ 


d iV xP 


dt 


(2.63) 


d^H 


d^M 


dP 


V^H + V[V-M) =£oiio-^ + £oiio-^-Vx 


(2.64) 


Onee more taking c = ——= and forming the wave operator on the left, 

a/^OMo 




1 d^H 1 d^M ^ dP ^^ 

^-v( v-M 

dp P dp dt 


(2.65) 


We have now obtained the wave equations in matter for E and H. 


^ d^E d^P / dM\ 1^,^ 
V^E-—^r^ =lio^r^ + lio Vx —-V V-F 


:2 dt^ dt^ 


dt y Co 


( 2 . 66 ) 
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(2.67) 


W^H- 


1 d^H 
dt^ 


1 d^M dP 

-Vx^r- VV-M 


dP 


dt 


This set of coupled hyperbolic inhomogeneous partial differential equations will prove 
resistant to the common methods of solution. Following an overview of transformation 
optics, which may be viewed as a method for circumventing a direct solution of the wave 
equations above for a desired propagation path, we will return to these expressions in an 
attempt to discover an iterative solution. 
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CHAPTER 3: 

An Overview of Transformation Optics 


We now proceed to a description of transformation optics, the technique currently used in 
the design of metamaterial structures capable of redirecting electromagnetic radiation as 
desired through those structures. 

As introduced by Pendry, Schurig, and Smith [2], transformation optics employs the in¬ 
variance of the Maxwell equations under coordinate transformations to convert the free- 
space wave solutions in a coordinate-transformed system, referred to as "electromagnetic 
space" [9], to equally valid solutions in physical space with electric permittivity and mag¬ 
netic permeability as determined solely by the coordinate transformation under considera¬ 
tion. A full understanding of the elegance, as well as the limitations, of the transformation 
optics approach will require a short overview of coordinate transformations as well as the 
covariant formulation of classical electrodynamics. 


3.1 Coordinate Transformations 

3.1.1 The Einstein Summation Convention 

Throughout this chapter, repeated indices indicate summation over all components up to 
the dimension of the space in which we are working. That is. 


n 

a%i=Y,a% (3.1) 

i=l 

where n is the dimension of the space, usually 3 (for classical mechanics) or 4 (for special 
relativistic applications). Furthermore, in order to avoid confusion between contravariant 
tensor components and exponents, should the need arise to express a tensor component 
raised to a power, the exponent will appear outside of a set of parentheses. 
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3.1.2 The Jacobian 


Consider an invertible coordinate transformation in four dimensions from some coordinate 
system into another coordinate system given by 


;c = ;c 


= x^ {x^,x^,x^,x^) 


(3.2) 


and similar for , x^\x^'. Coordinate differentials in the new coordinate system are given 
by 


, 3x^' 

dx^ = -^dx^ 
dx^ 


(3.3) 


The expressions form the elements of the Jacobian matrix 


““ dx^ 


(3.4) 


=A\dx^ 


(3.5) 


A contravariant vector, with components denoted obeys the above transformation law. 


y“'=A“«y“ (3.6) 

Just as the coordinate differentials comprise the primordial contravariant vector, the gradi¬ 
ent serves as the prototypical covariant vector. Note that for some function f. 


dx^' dx^dx^' 


(3.7) 


A covariant vector, with components denoted Va will obey the above transformation law. 
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3.1.3 The Metric Tensor 


The covariant metric tensor gjj is the collection of inner products of basis vectors in a given 
coordinate system. That is, 


gij = ti-tj (3.8) 

Note that due to the commutative nature of the inner product, gij is symmetric. In addition, 

ds^ = gijdx^dx^ (3.9) 

A metric gij therefore defines a geometry on the coordinate system in which we are inter¬ 
ested [16]. Covariant vector elements may be related to their corresponding contravariant 
vector elements through the metric: For some vector V, 


Vi = gijVJ (3.10) 

r = g'% (3.11) 

These operations are commonly referred to as “lowering” or “raising” an index, respec¬ 
tively. We may apply these relations to extract a useful relationship between the covariant 
and contravariant metric. 


Vi = gijVj = gijgj% (3.12) 

StVk = gijgj% (3.13) 


Where is the Kroneker delta. 
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)Vi = 0 


(3.14) 


Since V is arbitrary, 


SijS 


jk 




(3.15) 


and we eonelude that is inverse to ■ 

3.1.4 Tensors 

A eontravariant tensor of rank 2 is a geometrie objeet, denoted whieh transforms under 
the eoordinate transformation given by the Jaeobian in the following manner 

(3.16) 

while a eovariant tensor of rank 2 TcejS transforms as 

(3.17) 

This definition may be extended to eontravariant and eovariant tensors of higher rank by 
repeated applieations of the Jaeobian. For a eontravariant tensor of rank n, 


'ra'^a'2...a'n _ 4 j^2 j^n n-'aia2...a„ 
G- J ^2 * * * 


(3.18) 


while for a eovariant tensor of rank n, 


T, , , =a"',a“2 

a'j 0.2 " ^ 


GiG2-‘-Gfi 


(3.19) 


A eontravariant pseudotensor of weight w and rank n is an objeet with eomponents 
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which transform as [16] 


T” 


= |A| 


w ,al' ial' A% 
^ a2'"^ 


2-0x0 


2-..Cln 


(3.20) 


where |A| is the determinant of the Jacobian. Pseudotensors of weight ±1 are referred to 
as “tensor densities.” 


3.2 Covariant Electrodynamics 

3.2.1 The Field Tensors 

The covariant electromagnetic field tensor in Minkowski spacetime may be written 


/ 


^aj3 ~ 


Q As rZ 


B 


-Bz 

0 

Br 


f \ 

By 


-B, 

0 


(3.21) 


in the basis (ct,.r,y,z) with Minkowski metric 


Sap 


/l 0 0 0 \ 

0-100 
0 0-10 
\^o 0 0 -ly 


(3.22) 


The contravariant electromagnetic field tensor may be found through 


/ 0 




C 

0 


^ By 


^ -Ba, 


c 

-B, 

0 

Bx 


.iz\ 


B 


(3.23) 
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The contravariant magnetization-polarization tensor density (weight -1) is 



( 0 P,c 

PyC 

Pzc\ 

5Jt«/3 _ 

O 

1 

-Mz 

My 


-PyC M, 

0 

-Mx 


\-P,C -My 

Mx 

0 / 

The contravariant displacement tensor density (weight - 

1)2)«/3 [ 


o 

1 

— DyC 

-dA 

2)«/3 _ 

DxC 0 

-H, 

Hy 


DyC Hz 

0 

-Hx 


yD^C ~Hy 

Hx 

0 / 


(3.24) 


(3.25) 


The electromagnetic field tensor, magnetization-polarization tensor, and displacement ten¬ 
sor are related through 


2 ) a/3 


J_^a/3 _ 
Mo 


(3.26) 


3.2.2 The Maxwell Equations 

The free current density (weight -1) four-vector is 

Jj = (cp/,//) (3.27) 

From Gauss’ Law (2.27) in matter, 

Vb = pf (3.28) 
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diDx + d2Dy + — pf 


(3.29) 


^i-2)®i+<92-2)°^ + <93-2)°3 = -7? (3.30) 

c c c c ^ 


+ < 922 )°^ + = Jj 


(3.31) 


Since = 0, 


<9oD°® + ^12)01 ^ ^ ^^2)03 = yo 


(3.32) 


While from the Ampere-Maxwell Law (2.30) in matter, 


VxH 



(3.33) 


The x-component of the above relation is 


^2^1 d^Hy — Jf,x 3 " cdQDx 


(3.34) 


d2H^ - dsHy = 7} - ^o2)°^ (3.35) 

<922)21 + (932)31 = 7} - <9o2)iii (3.36) 

But 2)11 _ gQ 

(9i2) 11+ (922)21+ <932)31 =7}-(902)111 (3.37) 
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<9o2)°^ + ^i2)“ +<922)2^+<932)^^ =j} 


(3.38) 


Similarly, the y-component of the Ampere-Maxwell Law in matter may be expressed 


< 9 o 2)®2 + ^ 12)12 + ^ 22)22 + ^ 32)32 _ j 2 (3 39 ^ 

While from the z-component, 


<9o2)03 + ^12)^3 ^ ^^2)23 + ^32)33 = jj (3.40) 

(3.32), (3.38), (3.39), and (3.40) may be summarized, in the Einstein summation conven¬ 
tion. 


(3.41) 

We note that the above expression contains both Gauss’ law and the Ampere-Maxwell law 
in linear media. Turning to the two sourceless Maxwell equations. 


¥•5 = 0 


V X E -1- 


dB 


0 


Note that (3.42) may be expressed 


(3.42) 

(3.43) 


diBi -|- (9252 -|- d^Bj, = 0 


(3.44) 


<9i 523 -h d2F31 -h d3F]2 = 0 


(3.45) 
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While due to the antisymmetric nature of the electromagnetic field tensor (3.21), 


dlF32 + d2Fj3 + ( 93^27 = 0 


(3.46) 


From the x-component of Faraday’s law, 


V 


X £ + 


dB 


0 


(3.47) 


d2E'3 — d3E2 + cdQBi = 0 


(3.48) 


cd2Fo3 + cd3F20 + cdQF32 = 0 


(3.49) 


d 2 Fo 3 + <937^20 + <9 o 7^32 = 0 


(3.50) 


Once again we exploit the antisymmetric nature of the electromagnetic field tensor to obtain 


d2F30 + <937^02 + <9o7^23 = 0 


(3.51) 


Similarly, from the y-component of Faraday’s law we obtain 


< 9 i 7^30 + d^F 01 + dQF 1 3 = 0 


(3.52) 


diFo3 + d^Fjo + dQF3i = 0 


(3.53) 


And from the z-component of Faraday’s law 
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diFo2 + d2Fio + (9o^2i = 0 


(3.54) 


diF2o + d2Foi + dQFi2 = 0 (3.55) 

Equations (3.45), (3.46), and (3.50)-(3.55) contain all possible three-element cyclic per¬ 
mutations of the four spacetime indices, and may be encapsulated as 

^aFpa + d^Faa + ^aF afi = 0 (3.56) 


Or equivalently. 


da 


)^o 


(3.57) 

where is the four dimensional Levi-Civita symbol. The quantity -e^^^^Fyg = 

is denoted as the dual tensor to the electromagnetic field tensor. 


= 

0 


(3.58) 

We note that the dual electromagnetic tensor 





1 ^ Bx 

By 

Bz \ 


^a/3 ^ 

-Bx 0 

—B — 

c 

c 

0 

h 

C 

-ix 

C 

(3.59) 


1 

Co 

1 

F 

c 

0 ) 



1 ^ 1 


may be obtained from the components of through B ^ E, ^ B. Equation (3.58) 
contains both Earaday’s law and V ■ 5=0. Maxwell’s equations in matter may thus be 
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expressed in eovariant form as 


= Jj (3.60) 

(3.61) 


3.2.3 The Covariant Constitutive Relations 

The eonstitutive relations linking the primary fields to the material response fields are taken 
as 


P = eoXeE (3.62) 

M = XmH (3.63) 

for homogenous, isotropie, linear media, where the sealar values Xe and Xm are the eleetrie 
and magnetie suseeptibilities of the material under eonsideration. For anisotropie linear 
media, we take these suseeptibilities as seeond rank tensors, leading to the relations 

P = eQXe*E (3.64) 

M = Xm*H (3.65) 

where * indieates tensor eontraction. By eomponents, 

F = eoxi^Ej (3.66) 
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(3.67) 


The covariant representation of these relations may be given as [4], [10] 

^ (3-68) 

where are the components of the rank four contravariant susceptibility. We may 

conclude from the known nature of 971 as a tensor density of weight -1 and the tensor 
quotient theorem [16], that x must also transform as a tensor density of weight -1. The 
traditional permittivity and permeability of the material in question may be extracted from 
the susceptibility tensor density: Inserting (3.26) into (3.68) [10], 

2)“/^ = (3.69) 

Mo Mo 

Since is antisymmetric, it may be expressed [10] 

f = i (f“^ - F^“) (3.70) 

= l (3.71) 

Inserting this relationship into (3.69), 

®“'* = ^ -;t“'’'"')fpv (3.72) 

Defining as [10] 
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(3.73) 


^aPyS ^ J_ f^ay^pS _ ^Pr^aS _ ^apyd 
2/io V 

Note that the displacement tensor density elements are, in terms of this new susceptibility 
tensor, 


^^P =X^P^^F^y (3.74) 

We will take the above relationship as our covariant constitutive relationship. Note again 
from the quotient theorem that since D transforms as a tensor density of weight -1 that X 
must also transform as a tensor density of weight -1. 

We pause for a moment to observe that our constitutive relationship above is completely 
linear and therefore does not explicitly account for observed nonlinear effects e.g., self- 
focusing, frequency doubling. 

We may extract the traditional bianisotropic coupling coefficients, such as the permittivity 
and permeability from the elements of X. From, for example, the (ctjS) = (01) term [10], 

2)01 _ (3.75) 

= (3.76) 


+X^^^^Fjo + X^^'^^FnFX^^^^Fi2+X^^^^Fi3 

I y0120p _i_ v0121 p _l_v0122p _l_v0123p 

+ \ t20+-^ r21+-^ r22+-^ c 23 

+ X^^^^F30+X^^^^F31 +X^^^^F32+X^^^^F33 (3.77) 
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_^D +X^^^^B, 

C C C C z y 

+X0i20:^+x0i2i52-x®i235;, + x®i30^-x0i3i5,. + ;t0i32^^ (3.78) 
c c 


cD, = (XOIOI -x0“0) ^ + (X0102 _x0120) ^ + (;^0103 _;^0130) ^ 

+ (x0123_^0132^^^^ ('^0112_^0121^^^ 


Comparing this relationship to the common bianisotropic constitutive relationship [15] 


b = eE + 0 (3.80) 

B = hH+CE (3.81) 

where e, (^, /i, ^ are the spatial tensors which couple the electric displacement and mag¬ 
netic fields to the electric and magnetic fields, 

^ (3.82) 

Note that since the electromagnetic field tensor F is antisymmetric, our constitutive rela¬ 
tionship (3.74) requires that X“^^''be antisymmetric with respect to the pairs /iV. There¬ 
fore 


gif = (3.83) 

Similar calculations yield expressions for the remaining components of the permittivity 

[ 10 ] 
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(3.84) 


= 



3.3 Transformation Optics 

Having constructed the contravariant rank 4 modified susceptibility tensor density 
[15], and analyzed the transformation properties of tensors [16], we may now investigate 
the transformation properties of the material parameters contained within the susceptibility 
tensor. As a contravariant tensor density of rank 4 and weight -1, under a coordinate trans¬ 
formation described by the Jacobian the contravariant components of the redefined 
susceptibility X transform as [10] 

^a'li'ys' ^ (3.85) 

We now investigate the way in which the Maxwell equations in media transform under such 
a coordinate transformation. From (3.60) and (3.61), 

(3.86) 

= 0 (3.87) 

Applying the coordinate transformation given by the Jacobian with elements to the left 
hand side of the Gauss-Ampere-Maxwell law (3.84), 


= |Ar‘A" 


(3.88) 


But 
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^ dxl^' dx^ dx^^ 


(3.89) 


and so 




(3.90) 


Since the four-eurrent density is a tensor density of rank 1, weight -1, it will transforms as 


7 / = \A\-^A-',Jf 


(3.91) 


The Gauss-Ampere-Maxwell law therefore transforms as 


=JJ 


(3.92) 


and is thus form invariant under eoordinate transformations. Turning to the Faraday-Gauss 
law, 




(3.93) 


All of the Maxwell equations are therefore form invariant under eoordinate transforma¬ 
tions, as we might expeet of sueh a set of tensor relationships. A valid solution set 
^a/3^ 2)“^, to Maxwell’s equations in one eoordinate system will provide a valid 

solution set ^ in another eoordinate system aeeording to our transfor¬ 

mation properties [10] 




(3.94) 
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(3.95) 


^a'p'ys' ^ (3.96) 

The form invariance of the Maxwell equations and the resulting transformation properties 
of the modified susceptibility tensor density X is the concept at the heart of Transformation 
Optics. 


3.3.1 Spatial Transformations 

Turning now to the special case of transformations which do not involve the time: 

t' = t x=x{x,y,z) y=y'ixA,z) z=zixA,z) (3.97) 

Note that under this set of restrictions 

A^o = 1 A®' = 0 (3.98) 

for all i from 1 to 3 (spatial index). From (3.84) the components of the electrical permittivity 
transform as 




(3.99) 


From (3.98), the only a or 7 which survive the summation are a = 0 and 7 = 0 . 


y'/ = A \A\-U^'pA^'gX^I^^^ 


(3.100) 


But from (3.84), Reindexing, [10] 
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(3.101) 


3.3.2 Application: Cylindrical Cloak 

As our first investigation into the design of an electromagnetic redirection structure, we fol¬ 
low the development of a right circular cylindrical device in [9]. Let unprimed coordinates 
represent physical space and primed coordinates represent our electromagnetic space. We 
consider a coordinate transformation of the form 

r=- — -r' + a (3.102) 

b 

<)) = f (3.103) 

z = (3.104) 

ct = ct' (3.105) 

in cylindrical coordinates, where a and b are positive real numbers, with b > a. Such a 
transform maps the line / = 0 onto the circular cylinder r = a and the circular cylinder r' = 
b onto the circular cylinder r = b. Figure 3.1 illustrates such a coordinate transformation. 

Let us assume vacuum in the electromagnetic space. Then e''-’' = Eodj, and = fiodj,. 
As the transformation (3.102) does not involve the time, we may apply (3.101) in order to 
determine the transformed permittivity in physical space 

= \A\-^A\,A^ (3.106) 

The Elements of the Jacobian are 
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Figure 3.1: Transformation for Cylindrical Cloaking Structure. Source: [9] 


4“ - 


(3.107) 


And so 


40 

A Q, 


d {ct) 


(3.108) 



d{ct) 

dx‘' 


= 0 


(3.109) 


A' 


0 ' 


dx^ 


0 


(3.110) 


for all spatial indices i,i' = 1,2,3. Computing the remaining elements of this Jacobian, 


. dx (9(rcos(0)) 
dx' dx' 


(3.111) 
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(3.112) 


a1 


V 


d {rcos{^)) dr' <9 (rcos (0)) <9^' d {rcos{^)) dz' 
dr' dx' dcj)' dx' dz' dx' 


Ai 


V 


dr' 


b — a 
b 


r'cos +acos 


dr' d {rcos{^')) d^' 
dx' d^' dx' 


Let 9? 


b — a 
b 


A\, = 9?cos ^ — rsin 

C7yi 


dx' 


d ( V + y'^ 

A Y = Rcos [(j)') —- - — rsin (^') 


(9 (tan ^ 
dx' 


A^j, =Rcos ^ + rsin 


x'Hi + K 


A\, =Rcos^ sin^ 


1 _ dx^ _ (9 (rcos (<^)) 
dx^' dy' 


A 


1 

2 ' 


^ rr, / 

r— L/?r cos 
dy' 


(^') +acos (^')] 


A ^2' = 


(9 (cos (^')) 
dy' 


+ 9?cos (^') 


(9y 


(3.113) 

(3.114) 

(3.115) 

(3.116) 

(3.117) 

(3.118) 

(3.119) 

(3.120) 
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A 


1 

2 ' 


{Rr + a) 



+ 7?cos (0^) 



J d { ^/x'^+yA d ( A/y2+y2 


/y 


AV = -(/;/ + a)"^+^cos(f)"- 


A'y = - (;;/ + «) +/;eos (f) sin (f) 


A * 2 ' = — ^ cos (^') sin {ij>^) — cos {ij>^) sin (0') 


A^2' = (^ ~ ^) i^') (^0 


A 


1 

3 ' 


dx^ 


dx 


= 0 





{Rr' + a) 



dx' 


+ 7? sin (^') 


(9/ 

(9y 


(3.121) 

(3.122) 

(3.123) 

(3.124) 

(3.125) 

(3.126) 

(3.127) 

(3.128) 

(3.129) 

(3.130) 
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A^j, = - [Rr' + a] +7?sin (f) (3.131) 

= (- 

= ( 3 . 133 ) 

^2^^ ^ ^os (^') sin (^') (3.134) 


/„ / \ sin(0') „ . 

(7?r + a) —+ R sin 


(3.132) 


= ^7?—cos (^') sin (^') (3.135) 

-4". = [('"'+“) sin (f)] (3.136) 

( 3 . 137 ) 

Hv) d/ 

A\ = {Rr' + a)^^+Rsm{(^')— (3.138) 

A\, = (Rr' + a) + (3.139) 

A^, = [R/ + a) + +^sin(^')^ (3.140) 
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(3.141) 


(3.142) 

A^2' = J?sin^ (^') + -^cos^ (^') 

(3.143) 


(3.144) 

The remaining components of the Jacobian are somewhat simpler. For i' 

= 1 , 2 , we obtain 

a3 - -0 

(3.145) 

While 


a3 -1 

^ 3' - az' ■ ^ 

(3.146) 


Our Jacobian therefore takes the form 


/j?cos^(^') + ^sin^(^') (J?-cos(^')sin(^') 0 ^ 

A= (j?—^)cos(^')sin(0') J?sin^ (^') + ^cos^ (0') 0 

Vo 0 1/ 

The determinant |A| of this Jacobian is 


(3.147) 
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|A| = + ^sin^ (^O) (^0 "*" (^O) ~ 

(3.148) 


|A| = + cos^ (^') sin^ (^0 +sin^ (^')) “ (^0 (^0 


(3.149) 


|A| = ^ (cos^ (^') +sin^ (^0) +2^cos^ (^') sin^ (^') (3.150) 


1^1 = ^(cos^(<^’')+sin2(<^’'))^ 


(3.151) 


lAH^ 

r 


(3.152) 


Applying our transformation relationship (3.106) for purely spatial transformations, with 


e''^' = eo5;:, 


= \A\ ^A\,A^yeQd], 


(3.153) 


e^i = \A\ ^A^-,A^.,eQ 


(3.154) 


e'2' = eo|AriA'.,(A^)/ 


(3.155) 


e'2' = eo|Ari(AA^)'^ 


(3.156) 
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( 3 . 157 ) 


e = eol^l 


/ 


(^') + psin^ (^') ( 7 ?—^)cos(^')sin(^') 0^ 

( 7 ? —p)cos(0')sin(0') 7?sin^ (^') + pcos^ (^') 0 

0 0 1 / 

^Rcos^ {(j)') + y sin^ {(j)') [R — cos {(j)') sin {(j)') 
[R - fr) cos (^') sin (^') 7?sin^ (^') + f, cos^ (^') 
0 0 


X 


V 


e“=£o^((^cos2(f) + ^sin2 




cos' 


{(j)') si 


sin 



= £ 0 ^ U^cos^ (<))') + ^ sin^ (<)>') + ( 7?2 + ^ ) cos^ (<)>') si 


r 

^2 


sin 



.11 


= Co^ ( 7?^cos^ (cos^ +sin^ (^0) ■*" (^0 (^0 +cos 


£ 


11 




(3.158) 


(3.159) 


(3.160) 



(3.161) 

(3.162) 
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£12 = £0^ (r (r - I) cos 3 (<^') sin ( 0 ') + ^( r -^'^ cos (<^)') Sin 3 (<^') 

+ ^ (^ - 7 ) cos sin^ {(!>') + ^ (^ - 7 ) cos^ (<^') sin (<^>') ) (3.163) 

£12 = £0^ (t? (7? - cos ( 0 ') sin ((^') + ^ (^ 7 ? - cos {(j)') sin (<^')) (cos 2 (<^') + sin 2 (<^')) 

(3.164) 

£12 = £ 0 ^( 7 ?+^) (7?-^)cos(f)sin(f) (3.165) 

£12 = £ 0 ^ (^7?2 - cos {(j)') sin {(j)') (3.166) 

£i^ = 0 (3.167) 

£21 = £0^ (r (r - cos^ {(!>') sin {(j)') + ^( r -^'^ cos {(j)') sin^ {(j)') + 

r(^R — cos sin^ ^ ~ cos^ sin (^') ^ (3.168) 

£21 =£0^ 7) 7 7) +sin2 (^')) 

(3.169) 

£21 = £ 0 ^( 7 ?+^) (7?-^)cos(f)sin(f) (3.170) 
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+'^^sin"^ {^') + ^^08"^ (^0^ (3.173) 


^7?^sin^ [^') (cos^ +sin^ i^')) + {^') 

(()))+COS^ (<))') )^ 

(3.174) 


(3.175) 

g31 = £32 ^ Q 

(3.176) 

g33 _ ^ 

(3.177) 


Collecting the elements of e, 

/7?^cos^(^') + ^sin^(^') ^7?^ —cos(^')sin(^') 0^ 

e = eo^ cos (())') sin(())') T^^sin^ (<))') +^cos^ (<)>') 0 (3.178) 

Vo 0 1/ 
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in accordance with results claimed in [9]. Note that since r' —)■ 0 as r —)■ a, the susceptibil¬ 
ity components and 'x}'^ are expected to diverge. This divergence indicates a 

possible issue which may arise when performing numerical computations in order to deter¬ 
mine the response fields (see Chapter 5). Specifically, we anticipate the possibility that the 
value attained from those calculations may vary considerably depending on the fineness of 
our mesh. 

Converting to cylindrical coordinates in the basis (r, <j),z), 

£'■'/ = \A\^^ (3.179) 

= |A r ^ AEcartesian (A)^ (3.180) 

where A is the Jacobian from Cartesian to cylindrical coordinates, given by the transforma¬ 
tion 


r = V ^ ^ ^ 


-1 fy 


z — z 


(3.181) 


Computing the Jacobian, 



dj_ 

dxi 


^A = 


( ^ 

dx 

dx 
\ dx 


dr dr \ 
dy dz 
d^ dtj) 
dy dz 
dz ^ 
dy dz j 


(3.182) 


/ ^ 
r 

> 

r 

0) 


/cos(^) 

sin(^) 

o\ 


2L 



— sm{(j)) 

cos((/)) 

0 

jr2 

yl 


r 

r 

VO 

0 

l) 


\ 0 

0 

1/ 


(3.183) 


The permittivity in cylindrical coordinates is thus 
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i'f 



/cos(^) sin(^) 0 \ 

sin((/>) cos(<j)) Q 

r r 

Vo 0 1/ 


'^R2cos2(^) + ^sin^(^) 

(r^- 

J 

1 cos (^) sin (^) 

0^ 

(^^-^)cos(()))sin(<))) 

R^sin 

^(<^>) + ^COs2(())) 

0 

V 0 



0 

1 / 

/cos(<))) 

r 

o\ 



X sin(0) 

cos(0) 

r 

0 

(3.184) 


V 0 

0 

V 




/' i' 


/ 


^R^cos{^) R^sm{^) 0 \ 




r'" 

0 


s(0) 


0 


0 

1 / 



sin(0) 

r 

cos((j)) 

r 

0 


o\ 

0 

1 / 


(3.185) 


e^j' I = 


r' 


fR^ 0 0 \ 
0 0 
VO 0 i/ 


[r'R 0 0\ 


= Co 


0 pV 


Vo 


0 


0 


(3.186) 



0 


eo 


0 0 


M 

0 


r—a / 
/ 


(3.187) 


Note that these results are not in accordance with [9], although both results are diagonal in 
cylindrical coordinates. We attribute the discrepancy (of a factor r) to the author’s neglect 
of the tensor density nature of the permittivity and subsequent omission of the necessary 
factor = r. Regardless, the permittivity is diagonal in cylindrical coordinates, as we 
might have predicted based upon symmetry requirements. Again we see that at least one 
component of the susceptibility is expected to diverge at r = a. 
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CHAPTER 4: 

An Iterative Approach to the Wave Equations 


In this chapter we will discuss an alternative approach to the wave equations in matter 
as determined in chapter two. We will then consider this approach when applied to the 
problem of scattering of a pure plane electromagnetic plane wave from a structure of known 
geometry and electromagnetic susceptibilities. 
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frequency domain. 


v^E + k^E 


-—P-icouo (v x^'] - —V (v-P 
So V / Co V , 


(4.5) 


where k = — and we have used the relationship c =-. 

c VeoMo 

frequency Fourier transform of (4.2), 


Similarly, from the time—)■ 


co^ 


co^ 


+ = —^M + iO)[WxP 


V V-M 


(4.6) 


+ k^ti = -k^§ + i(o {y Xpj-W {w 


(4.7) 


where k is again equal to —. We recognize (3.5) and (3.7) as a set of inhomogeneous, 

c 

coupled Helmholtz equations for the components of E and El. 


4.2 Constitutive Relations 

We wish now to introduce a set of constitutive relations between the polarization and mag¬ 
netization fields P and M, and the electromagnetic fields E and H, respectively. In linear, 
isotropic media, the relations 


P = eox'E 

where x^^ x'^ the electric and magnetic susceptibilities, and are employed to de¬ 

scribe the response of the polarization and magnetization to the fields. In order to account 
for the possible anisotropy of our media, we express each of the susceptibilities as the 
second-rank tensors X^ X^- The constitutive relations now read 
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(4.8) 


P = eoX‘'*E 

^ = (4.9) 

where * denotes tensor contraction. By components, 

fa = £« I X’afEj, (4.10) 

(i.l 

«a=i4“o«;j (4.11) 

/3=1 

Our Fourier-transformed wave equations now read 

V2| + = -k^ *£)- i(OHo (Vx -V(^V- (^x" *^)) (4.12) 

+ k^^ = -k^ l^x "’*+ io}eo (Vx (^x"*^))-y(y- (x "’*) (4.13) 

4.3 Incident and Response Fields 

We now consider separately the incident fields and the material responses to those fields. 
That is, 

P=£o + Pr (4.14) 

n=nQ+nr ( 4 . 15 ) 
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where Eq, Hq are the incident fields and Ey, Hy the fields due to the material response. From 
(4.12), 

- * (io + i) ) - ( V X [x"' * (^4 + - V (v ■ * (4 + i) ) ) 

(4.16) 

Let us model the incident fields as arriving through vacuum. Then the incident fields are 
expected to obey the sourceless Helmholtz equations: 

v2fo + A:2fo = 0 (4.17) 

Thus, 

V^iy + k^iy = 

- * (^0 + i) ) - ( V X (^x" * (4 + )) - V (v ■ * (4 + i) )) 

(4.18) 

Similarly, from (4.13) 

W^ffy + k^ffy = 

- e [x'^ * (4 + ^r) ) + icoeo (V X [x^ * (^4 +_ V (v ■ [x"' * (4 + A) ) ) 

(4.19) 

Note the near symmetry in the source terms of the above two equations for the response 
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fields. We define 


Fi{x,Y^=-k^{x‘'*^)-i(OliQ{yy< {x"'*y))-y{^-{x^*X^) (4.20) 

F 2 (^, ?) = -k^ [x" *+ i(oeo (v x * ?) ) - V (v ■ [x'^ * i) ) (4.21) 

Note that 

Fi (Xi+X2ji+Y2) = 

- * (^1 + ^2) ) - i(OFo (v X (^x" * (i^i + >^ 2 ) ) ) - V (v ■ [x^^ * (^1 + ^2) ) ) 

(4.22) 

= -k^ (z" *^i) - (x‘'*X2) - icoFo (v X * fi) ) - icono (v x (^x" * Y2) ) 

-V (v- (;t^*Xi)) - v(v- [x‘^*X2)) 

Fi (Xi+X 2 j,+Y 2 '^ =Fi (Xi,yi) +Fi (^ 2 , 72 ) (4.23) 

And similar for F 2 . From (4.18) and (4.19), 

v2|, + k^i, = Fi (4,4) + (4, 4) (4.24) 

¥^4 + ^24 = F2 (4,4) + i^2 (4,4) (4.25) 
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4.4 An Integral Equation 


The Green Function for the Helmholtz equation in two dimensions is known (see Appendix 
A) to be given by 


Gyto (g r') = (ko I r - /1) (4.26) 

9 (1) 

For the Helmholtz operator + Uq, where Hq is the first Hankel function of the first 
kind. With this Green function in hand, we may express the solution to our inhomogeneous 
Helmholtz equations (4.24) and (4.25) as 


i = Gk (F, r') (Fi (fo (r ), 4 (r)) + F, (i (r') , (?')) ) {All) 

Gk (F,F') (f 2 (4 (F') Jo{r))+ F 2 (A (F'),i (F'))) d^r' (4.28) 

We have now converted our differential equations for the response fields in equations 
(4.24) and (4.25) into a set of coupled integral equations for these fields. In an attempt to 
keep our expressions under control, we further introduce the notation 


F[{xj) =Fi i(F),y(F) 


(4.29) 


and similar for F 2 . Our integral equations for the response fields now read 


ir= / Gfc(F,F')Ff (^0,4) / G4F,F')Ff (irA) d 


2 / 


(4.30) 


/ Gk (F, r') (A ^ 0 ) d^r' + / G, (F, r') f/ A^iAd 


2p/ 


(4.31) 
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4.5 The Neumann Series 


We propose to attain an approximate solution to equations (4.24) and (4.25) for the response 
fields through the following iteration method: Insert the entirety of equation (4.27) and 
(4.28) for Er and Hr, respectively, back into the last terms of (4.27) and (4.28). This method 
produces what is known in mathematics as a Neumann series, which is often employed in 
the approximate solution of integral equations [17]. The technique is also commonly used 
for quantum scattering problems, where truncation of all but the first term of the resulting 
series is known as the Born approximation [18]. 


72-/ 


Er = J_^Gk{7,r)F[ [^Eo,Ho) d 

+ rGk{r,r')Ff( ( T (44) + I Gk{r,r')F{' (44) d 


Gk{f,r)F[ (4,4) d^r+ I Gk{r,f')F^' (4,4) d 


;/ -//^ 177" 


72-// 
7 ^ 


72-// 


4r (4.32) 


4 = Gk (F, r') 4 (4,4) d^f 


/ oo 

4(F,F')4 

'CXD 

/ OO 

4(F,F')4 

'CXD 



(4,4) Gk{7',r)F{" (4,4) d^r^ d'^r' 

(4,4) G 4 r', 4 ) 4 " (4,4) d ^^ JV 


(4.33) 


where we have used the linearity of the function Fi from equation (4.23) in the last step. 
Provided the convergence of the resulting series, we may continue this iteration process 
as we wish in order to attain an approximate solution for 4- Similarly, for the response 
magnetic field 4, 
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^r = I^^Gk{7,r)F^ (4,4) d^r 





j2^// I 

a r + 




5 



4,4 4'(4,4) 


d^r+ 


Gk 


f'y 


)4 


F H 

lly 



(4.34) 


4=/ G4r,/)F|'(4,4)^ 

J —CX3 

/ CO 

G,(?.?')f2'' 

-CO 

poo 

+ Gk{f,f')F^ 


]2p 


Gk{f,r")Fr iHo,Eo) d^r, I Gk{f,r)Fr (Eo,Ho) JV' jV 


Gk (r , 7”) F|" (HnEr) d^r" , I Gk (r ,r') Ff" (F„j 2 ^/ 


(4.35) 


Identifying the first term of (4.33) and (4.35) as the first-order response in x, second term 
as the second-order response, and the final term the remaining higher-order response, we 
may continue the iteration process as we desire by inserting each of (4.33) and (4.35) into 
the expressions for the response fields of each of these expressions, each time attaining the 
next higher order solution to our wave equation. 


4.6 Incident Plane Wave 

As a precursor to our implementation of the iteration scheme outlined above (see Chapter 
5), let us consider the material response of a two-dimensional structure embedded in the x-y 
plane when subject to an incident electromagnetic plane wave, polarized in the z direction, 
with angular frequency ft)o, travelling in the x direction, and arriving through vacuum. 


Eq = 


(4.36) 
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where kQ is the wave veetor, oriented in the direetion of propagation, with magnitude given 

K 7 ^ 

by ko = —. 
c 


4.6.1 First Order Response 

The temporal Fourier transform EQ(r, co) with respeet to waves propagating outward is 

Eo= (4.37) 

J —oo 

EQ = EQe^^^ (4.38) 

J —oo 


So = iKEoe^^’^d {(0-(0o)z (4.39) 

where 5 (to — tOo) is the Dirae delta distribution. For plane eleetromagnetie waves, Hq is 
determined by our ehoiee of polarization and propagation direetions to be 




(4.40) 


Where//o 

wave is 


^0 

MoC V /io 


EocEq. The temporal Fourier transform of this ineident plane 


So = {a)-(Oo)y (4.41) 

We now take x‘^ = iii order to ensure that the strueture is impedanee matehed to free 
spaee, eliminating refleetion - a desired property in a CDEW applieation. The first-order 
response eleetrie field is thus, from (4.33), 


ErA = 


Gk{r,S)E{' (So,So) d 




(4.42) 
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Eni = 

x(^-k^(^X*eo (/)) - ico^o (Vx(^x* 4 {r) ))-v(v-(^x*eo{r)))) d^r (4.43) 

X (k^ (^Xe * 4 {r) ) + ico^o (v X (^Xm * 4 {r) ))+v(v-[xe*eo (/)) ) ) d^r' (4.44) 

Inserting the expressions (4.39) and (4.41) for £0 and Hq, respeetively, and exploiting the 
relation Xe = Xm for nonrefleetanee, 

l,,l =-5(®-tOo)y y H^^\k\r-f\) 

{x * z) ) - /ft) Mo (V X j) ) + V (v ■( z) ) ) ) d^V 

(4.45) 

Er,l^-5{(0-COo)^ J H^^\k\r-r'\) 

[k^Eo (^x * z)) - iCOI^oHo {^x{x* y)) + EqV {v ■ [x * z) ) )) d^r' 

(4.46) 


Let us express x in tho basis {x,y,z} as 
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IXn Zi 2 0 

X = X 21 X 21 0 

\0 0X33 


(4.47) 


since all of the TO-derived suseeptibilities under eonsideration in this work take the above 
form. Furthermore, sinee all of the eomponents Xij have no z-dependenee, the tensor eon- 
traetions of equation (4.46) are 


X* e‘^^^z)=e‘^^^X33Z 


(4.48) 


V-(x*e‘^^^z) =V-(e‘^^^X33z) = 




(4.49) 


iko^ A _ JkQX ^X33 ^ 




(4.50) 


As stated, ;t ^33 does not vary with z, and so 


V- (x*e'’^^^i)=0 


(4.51) 


Turning to the eurl term in (4.46), 


X * {X\2^ + X22y) 


(4.52) 


Tvy ( f i/tnxAM <5 1+;^:225/ 2 )) 

^x[X* ° y))\ . = £jki— - ^^ - — 


(4.53) 
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d (g'^°^Z22) A d ^ ( d _ d (g^'^°^Zi2) A . 

dz j dz I dx dy I 

(4.54) 


Once more exploiting the z-independence of the susceptibility components, 


V X (z * () = I 


(4.55) 


ifcoxx)) _ I Jkpx ^X22 ikQXiie'^'' - ^ z 

dx dy J 


Vx[x*[e^^^y)) = (e 


(4.56) 


Inserting (4.48), (4.51), and (4.56) into our expression for the first-order response electric 

E E IE 

field (4.46), taking 5o = — —)■ UqHq = — —> //q = \ —Eq, 

c c \ llQ 


Er,l = -5 (to - too) 


nEoi 

~Y~ 




. 0 ) 

X 33 - 1 — 

c 



+ ikoXii 



(4.57) 


But k 


(0 

—, so 
c 



-5 (to- too) 


nEoi 

~Y~ 




X 33 - ik 


dX22 

dx 


+ kokX22 + ik^^^^ zd^f 


(4.58) 


4.6.2 Second Order Response 

Returning to equation (4.33), the second order response electric field is 


42 =/ Gk{ry)E((l G,(/,/')< (4,4) 


Gk{r',r)EnHo,Eo)d^r]d 


(4.59) 
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The first integral over r is identical in structure to the first order response field, given in 
(4.58). 

Gk (/,/') Ff ( 4 , 4 ) = 

-d{co-coo) J 4^^ ~ I) (^^^X33 - + kokXii + 

(4.60) 

Focusing on the second integral over r', 

I^^Gk{r',r)Ff (4,4) = 

Gk (/, /') [-k^ (z * 4) + /®eo (v X [x^ * 4) ) - V (v ■ [x’" * 4) ) ) dV 

(4.61) 

= 2;r5 (to - o^) Gk (F',F") (^ - * (-^oj) + ^®eoV x {x^' * {Eqz)) 

-V(V.(;t'”*(-^oy))))t^V' 

= ^^5(to-G^) [ H^^\k\f-r"\)F’'°'^'(k^.[^x*y + iO}eo^xix'*z) 

L J —oo \ y jUq 

+ A/f^v(v-(r*y))V"^" 

V Mo / 
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nEoii 


S(co-coo) [ Hl^^\k\r-r"\)e''^^''(k^.f^{Xi2X + X22y) + icO£oyxiX 33 z) 
7—00 y V A^o 


+ J^^{^-{Xnk + X22y)) jdY 
V Mo 


TTEgi 


S (co- 


cOo) r {k\f-r\)e‘’^°^" 
J —oo 



— ixnx + X22y) + ifoeo 


dX33 . dX33 


dy 


dx 


-y 


(dxn dx22 
^Vmo <?/ 


(fr 


Gk{k',r)E^' {HQM(fr = 


kEqI 


5 {co 


-®o) / Hq'’ {k\r -f'\) (k^.f^ {Xi2X + X22y) + iO)£o ( 

7—00 \ y [1 q \ 


dX33 . < 5;^'33 


a/ 


-;c — 


dx" 


£o 

+ W — 

Mo 


^ <?^Z22 


a(/)2 <9/^/ 




<9^Zi2 ^ < 92^22 


[<9/<9/ d{y")\ 


y (fr (4.62) 


Inserting (4.60) and (4.62) into (4.59) yields 
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Gk (r,r) X 



PI 


kEqI 


5 (w - gjd) j {k \f -r"\) (^^X33 - + hkXii + ’ 


5 {co 


CDo) r {k \r' -r\) {k\l^ {Xn 

J—oo y V 


— (Zl2-^ + Z22y) + i®e0 


< 5^33 . dX33 . , , 
3C- zr-^y + 


dy 


dx 



d^Xn ^ d^Xii 


Po\[d{x) dx dy 


x + 


d^X 

dydx" • d{/y 


12 d^X22 


d^f IjV (4.63) 


Er.2 = 


Gk (r, r ) 5 (ft) - a»o) \k^ ' 


■X 


X* 


(^J H^'’ [k\r -r”\) (k^X33-ik^^+ kokX22 + ik^^^ zd^r" | | + 


nEo 


co/ioV X 


m 


^\k\r-r"\)e {Xnx + X22y)^ + 


/ft) Co 
kEqI 


f dX33 dX33 




+ 



d^Xn ^ d^X22 


dy'dx d{yy 


n + 


V (v ■ (;t * (^"^33 - ik 


., <5;^f22 


+kokX22 + ik^^^ 


zd r 


(4.64) 


We may observe from the relative complexity of the above result that the second and higher 
order terms of the response fields (4.33) and (4.35) will be computationally taxing for all 
but the simplest structures and material parameters. We will therefore focus our efforts in 
the following chapter on characterizing first order responses for the TO-derived structures 
which appear in the literature. 
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CHAPTER 5: 
Applications 


In this chapter we will apply the iterative approach outlined in Chapter 4 to some basic 
electromagnetic redirection structures and analyze the results in an attempt to quantify the 
error fields which result from nonlinear effects in these TO derived devices. 


5.1 Cylindrical Shield 

We will attempt to apply the iterative method of Chapter 4 to the cylindrical structure 
presented in [9] in Cartesian coordinates, followed by an analysis of the same structure in 
the coordinate system most suited to such a structure, cylindrical coordinates. 


5.1.1 Cartesian Coordinates 

From equation (3.178) the electric permittivity e is, in the basis {x,y,z). 


J 

( cos 

2(<))) + ^sin2(<),) 

(^^-^)cos(<)))sin(<))) 

0^ 

r 

£ = £o — 

Rr 

{r^- 

^)cos(<)))sin(<))) 

7?^sin^ {(j)) + ^cos^ {(j)) 

0 


V 

0 

0 

1/ 


(5.1) 


The electric susceptibility Xt corresponding to this permittivity is 


Xe 


1 

—£ - 

£o 


/l 0 0\ 
0 1 0 
Vo 0 1 / 


(5.2) 


Xt = 


'^^cos^(^) + ^sin^(^)- 1 

V 0 


^ sin^ (^) + ^ cos^ (^) - 1 
0 


0 ^ 
0 



(5.3) 
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(5.4) 


cos2 (<),) + ^ sin2 (<),)- 1 (^ - ^) cos {(l>) sin (<))) 0 

Ze= (^-^)cos(<^)sin(<^) ^sin2(())) + ^cos2(()))-l 

V 0 0 ^-1 

where we have used the transformation relation (3.102) in the last step. The permittivity is, 
as for the case of the square structure of [8], of the form 

fXn Zi2 0\ 

Xe= Xi\ Xii 0 (5.5) 

VO 0 X 33 J 

and so our result (4.58) for the first order response to an incident plane wave, polarized in 
the z-direction will apply to this cylindrical structure as well. 



Er,l = 5 (to — tt^) 


nEoi 

~ir 




X33 - ik 


dXll 

dx 


+ kokX22 + ik 


dXn 

dy 


(5.6) 


Constructing the permittivity as per [9] (see MATLAB code. Appendix B) for a structure of 
inner radius 1 meter and outer radius of 1.5 meter, we note immediately from Figure 5.1 the 
sharp peaks at the inner radius (set to 1 meter in this simulation). The mathematical form of 
the electric susceptibility from equation (5.16) suggests that these peaks are an artifact of 
the rectangular grid enforced upon an inherently cylindrical structure. That the location and 
height of these peaks, in addition, varies with mesh size (see Figure 5.1), provides further 
evidence of a granularity problem. This granularity issue and the resultant inablility 
of this model to accurately determine susceptibility values according to the TO 
prescription suggests that any application of the iterative approach outlined in Chapter 4 
which uses these susceptibility values should be suspect. As seen in Figure 5.2, the 
attained first order response field does not exhibit any of the expected behavior of a TO- 
derived redirection structure. We attribute this incongruence to the issue of representing a 
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a circular structure on a rectangular grid. 


for Cylindrical Kedirection Structure 



y (mrters) .2 .2 x (meteix) 


Xfj. for Cylindrical Redirection Structiire 



Figure 5.1: Left: xx component of the electric susceptibility as computed in cartesian coordinates, 
for a redirection structure of inner radius 1 meter, outer radius 1.5 meter, and mesh size of 0.05 
meter. Right: The same simulation, with a mesh size of 0.03 meter. Note the sharp peaks which 
arise due to the granularity of the simulation regardless of mesh size. 


1st Order Response Field Amplitude Relative to \E{)\ 


1st Order Combined Field Amplitude Relative to lifnl 



Figure 5.2: Left: First order response field for incident ^-polarized plane wave with wavenumber 
1 X in Cartesian coordinates, mesh size 0.03 meter Right: Combined response field. 

Note that the jagged peaks of Figure 5.1, an artifact of the Cartesian grid imposed on the 
structure, have carried over to the response field 


5.1.2 Cylindrical Coordinates 

In an effort to address the granularity issue whieh has likely eompromised the usefulness 
of the treatment seen in the previous subseetion, we turn to the eoordinate system naturally 
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suited to a cylindrical structure — cylindrical polar coordinates. Constructing a cylindrical 
coordinate grid and inputting the susceptibility values determined from (3.189) and [9] (see, 
once again, MATLAB code. Appendix B), we observe that the sharp peaks characteristic 
of the susceptibility in a Cartesian grid system appear to have been eliminated in favor of 
the expected smoothly varying values from (3.178), as we may see in Figure 5.3. 


for Cylindrical Uedirection Structure 

60 

40 

S 20 

X 

0 

-20 
2 

y (meters) .2 .2 x (meters) 

Figure 5.3: xx componenet of the electric susceptibility for the same cylindrical structure as seen 

2k 

in Figure 5.1, in cylindrical coordinates. Mesh size is 0.03 meter in the radial direction, - 

, 200 
radian azimuthally. 

From equation(4.47), for an incident plane wave polarized along the 2-direction, 



Er,l = S{C0-COo)^ f (k|r-r I) 

Z J — OO 

z) ) - iCO^oHo (v x y) ) + EqV (v-(^x* z) ) ) ) d^r 

(5.7) 

Performing the computation (5.9) with this cylindrical mesh, we note that although the 
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overall combined (incident and first order response) electric field amplitude shows overall 
promising behavior, falling rapidly within the structure, the resultant fields can be made to 
vary over several orders of magnitude depending upon the choice of mesh fineness. See 
Figures 5.4-5.6 for a comparison of the combined electric field for various mesh sizes. 


1st Order Combined Field Amplitude Helative to |£'|)| 



y (meters) -2 2 

X (meters) 



Figure 5.4: Left: Surface plot of the combined (incident and first order response) electric field 

2ti 

amplitude for the cylindrical structure. Meshsize is 0.003 meter in the radial direction, 

azimuthally. Wavenumber is 0.01 Right: The response field for the same meshsize, along 
the x-axis. 


ist Order Combined Field Amplitude Relative to |£’|)| 1st Order Combined Field Along X-axis Relative to \Eo\ 
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Figure 5.5: Surface plot and two dimensional plot along the x-axis for the same structure as in 
Figure 5.4, with a meshsize of 0.005 meter in the radial direction. 
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1st Order Combined Field Amplitude Relative to |£‘ii| 
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Figure 5.6: Surface plot and two dimensional plot along the x-axis for the same structure as in 
Figure 5.4, with a meshsize of 0.007 meter in the radial direction. 


In addition, our results appear strongly dependent upon our ehoiee of ineident wavenumber, 
as shown in Figures 5.7-5.10. 
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Figure 5.7: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the cylindrical 
redirection structure simulation, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber O.OOlm^F 
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1st Order Combined Field Amplitude Relative to |£‘ii| 



y (meters) .2 .2 x (meters) 



Figure 5.8: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the cylindrical 
redirection structure simulation, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber O.Olm^^ 
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Figure 5.9: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the cylindrical 
redirection structure simulation, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber 


We note from Figure 5.10 that as the wavelength falls to less than three orders of magni¬ 
tude relative to the radial mesh size, the simulation appears to break down altogether. We 
attribute this failure to the inability of the (now relatively eoarse) mesh to fully extraet the 
neeessary data from our model for the ineident eleetromagnetie wave. 
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Figure 5.10: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the cylindrical 
redirection structure simulation, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber . 


5.2 Reduced Material Parameter Cylindrical Shield 


In an effort to eliminate the influenee of diverging material parameters upon our results, 

we turn to an approaeh diseussed in [3] and [11]. Enforeing = 1, Cummer et al. 

£‘j a'j 

have shown that although the non-refleetanee eondition — = — is no longer met, the 

£o jUo 

remaining relevant material parameters in eylindrieal eoordinates reduee to 



(5.8) 



(5.9) 


The absenee of eleetrieal permittivity eomponents other than will prove a useful sim- 
plifieation when eonsidering the strueture’s response to a plane wave polarized along the z 


direetion. 
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We note that these redueed material parameters involve only one eylindrieal eoordinate 
component which varies with position through the structure, and that none of the com¬ 
ponents of the prescribed susceptibilities diverge anywhere within the structure (see Fig¬ 
ure 5.11). These reduced parameters are therefore expected to be advantageous from a 
manufacturing standpoint [11], as arbitrarily large electromagnetic material parameters are 
not currently realizable. In simulations, the reflected intensity arising from now nonzero 
impedance mismatch between surroundings and structure is somewhat minor [11]. 


Xxx for Cylindrical Redirection Structure 


Xyy for Cylindrical Redirection Structure 






Figure 5.11: Xm.xx and Xm,yy for the reduced material parameter cylindrical structure. Unlike 
the material parameter solutions of the previous section, these susceptibiity components do not 
diverge as we approach the inner boundary. 


Since, for this application, Xe 7 ^ Xm, the general first order response (4.44) becomes, fol¬ 
lowing development identical in form to that of equations (4.45) to (4.58), 



-d {(O-COq) 


nEoi 

~Y~ 




Ze,33 


+ kokXm,22 + 


dx 


dy 
(5.10) 


Implementation of this expression for the first order response electric field in MATLAB 
(Appendix B) and analysis of the resulting combined (incident and first order) electric field 
demonstrates that the previously observed large variation of results with respect to meshsize 
has been eliminated. See Figures 5.12-5.14 for a depiction of this simulation stability with 
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respect to meshsize. 
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Figure 5.12: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the MATLAB 
simulation of the cylindrical reduced material parameter redirection structure, inner radius 1 
meter, outer radius 1.5 meter, with cylindrical mesh, radial meshsize of 0.05 meter, azimuthal 
meshsize and wavenumber 


1st Order Combined Field Amplitude Relative to |£‘||| 


0.96 

0.94 ^ 



y (meters) .2 .2 x (meters) 



Figure 5.13: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the MATLAB 
simulation of the cylindrical reduced material parameter redirection structure, inner radius 1 
meter, outer radius 1.5 meter, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber O.lm^^ 
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Figure 5.14: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the MATLAB 
simulation of the cylindrical reduced material parameter redirection structure, inner radius 1 
meter, outer radius 1.5 meter, with cylindrical mesh, radial meshsize of 0.008 meter, azimuthal 
meshsize and wavenumber 


We still observe, however, a very strong dependenee of our results on the ehoiee of ineident 
wavenumber. See Figures 5.15-5.19 for a review of simulation dependenee upon ineident 
wavenumber. 
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Figure 5.15: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the MATLAB 
simulation of the cylindrical reduced material parameter redirection structure, inner radius 1 
meter, outer radius 1.5 meter, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber O.Olm^F 
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Figure 5.16: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the MATLAB 
simulation of the cylindrical reduced material parameter redirection structure, inner radius 1 
meter, outer radius 1.5 meter, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber 
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Figure 5.17: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the MATLAB 
simulation of the cylindrical reduced material parameter redirection structure, inner radius 1 
meter, outer radius 1.5 meter, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber . 


As previously, the simulation appears to breakdown as the ineident wavelength falls to 
within three orders of magnitude of the radial meshsize. 
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Figure 5.18: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the MATLAB 
simulation of the cylindrical reduced material parameter redirection structure, inner radius 1 
meter, outer radius 1.5 meter, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber lOm^^ 
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Figure 5.19: Surface plot and two dimensional plot along the x-axis of the amplitude of the 
combined (through first order) electric field, relative to incident field amplitude, for the MATLAB 
simulation of the cylindrical reduced material parameter redirection structure, inner radius 1 
meter, outer radius 1.5 meter, with cylindrical mesh, radial meshsize of 0.02 meter, azimuthal 
meshsize and wavenumber lOOm^^ 


Based upon the above results, we hypothesize that the strong dependence on meshsize ob¬ 
served in Section 5.1 is a result arising from the divergence of the material parameters 
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observed in those structures. As our mesh becomes finer, it incorporates a greater portion 
of the region near the inner boundary, where susceptibility values diverge, leading to rising 
response field values, as we have seen. The simulation results obtained in the case of the 
reduced parameter structure, which do not involve diverging electric or magnetic suscepti¬ 
bility, appear highly stable with respect to meshsize variation. A common feature to both 
approaches, however, is the strong dependence of our results upon incident wavenumber. 


5.3 Square Shield 

We now consider the square electromagnetic redirection structure outlined in [8]. Note that 
as for the unreduced parameter cylindrical structure, some components of the susceptibility 
diverge as we approach the inner surface. In Cartesian coordinates, the susceptibility tensor 
X may be expressed as 



(Xn 

X\i 

0 \ 


Xi\ 

Xii 

0 


^0 

0 

^ 33 / 


(5.11) 


For an incident plane wave polarized in the z-direction, the first order response is given by 
equation (4.58) as 


Ey^x = d((0-coo) 


kEqI 






k^X33 


+ kokxii + 


zd r 


(5.12) 


Results of simulations (Appendix B) for this structure are similar to those observed for 
the Cylindrical structure once an appropriate mesh shape was imposed upon it. We again 
observe similar overall behavior between all structures, displaying some of the desired 
combined field behavior (field amplitude falls to a minimum rapidly upon entering the 
inner area of the structure), with amplitudes strongly dependent upon meshsize, as shown 
in Figures 5.20-5.22. 
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1st Order Combined Field Amplitude Relative to |£‘ii| 




Figure 5.20: Surface plot and plot along the x-axis for of the combined (to first order) electric 
field amplitude, relative to incident field amplitude, for the square structure simulation, with 
meshsize 0.03 meter, wavenumber 0.02 half side lengths 1 and 2 meters. 
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Figure 5.21: Surface plot and plot along the x-axis for of the combined (to first order) electric 
field amplitude, relative to incident field amplitude, for the square structure simulation, with 
meshsize 0.05 meter, wavenumber 0.02 half side lengths 1 and 2 meters. 
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Figure 5.22: Surface plot and plot along the x-axis for of the combined (to first order) electric 
field amplitude, relative to incident field amplitude, for the square structure simulation, with 
meshsize 0.07 meter, wavenumber 0.02 nT^, half side lengths 1 and 2 meters. 



In addition to these findings, we observe eombined field amplitudes dependent on 
wavenumber as per the eylindrieal strueture. Also noted is the breakdown of the simu¬ 
lation similar to that of Figures 5.10, 5.18, and 5.19 upon the ineident wavelength lowering 
to within three orders of magnitude of the meshsize. 
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CHAPTER 6: 
Conclusions 


Here we will consider the results of the approach to electromagnetic scattering pursued 
through Chapters 4 and 5, and outline possibilities for future work. 


6.1 Effectiveness of the Iterative Approach 

Although well-grounded in the fundamental relations of classical electromagnetism, im¬ 
plementation of a Neumann series approximation for the response electric and magnetic 
fields can be problematic when applied to exotic structures such as TO-based CDEW elec¬ 
tromagnetic redirection structures. The divergence in the material susceptibility required 
for even the simple ideal structures we have examined here can lead to considerable vari¬ 
ation in the results of our numerical calculations depending on mesh size and incident 
wavenumber. For this reason, it is likely that, in order to attain convergence of numeri¬ 
cal results, a meshsize finer than possible on the hardware employed for this study should 
be used, or analyses should be restricted to reduced material parameter structures such as 
those studied in Section 5.2. Regardless, powerful hardware would also be required to be¬ 
gin analyzing the numerical results for second and higher order responses. These higher 
order results would be of value in addressing the remaining issues of series convergence 
and the extreme incident wavenumber dependence of our results. 

Due to practical engineering concerns, a real-world CDEW redirection structure will fea¬ 
ture neither a continuously varying susceptibility, nor a susceptibility with any components 
diverging as a particular bounding surface is approached. Studies such as [3] and [13] have 
examined the linear response of layered structures, designed to approximate the cylindrical 
redirection device proposed in [9] and analyzed in Chapter 3 above. 


6.2 Future Work 

Opportunities for follow-on research are plentiful. As mentioned, attaining a convergent 
first order response and analyzing the numerical computations for second and higher order 
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response fields may be possible with aeeess to greater eomputing power than was available 
for this work. In addition, onee the higher order responses have been well eharaeterized, 
and assuming that the Neumann series is shown to eonverge, the effieaey of the approaeh 
may be tested via experiment using the layered struetures whieh have been shown to elosely 
mimie the theorized linear behavior of TO-based deviees. Sueh experimental work eould be 
earried out as an extension to the metamaterials researeh eurrently being done by Professors 
James Luseombe and Dragoslav Grbovie of the NPS Physies Department. 

Beyond these refinements, there is eonsiderable theoretieal ground whieh may be analyzed 
through the approaeh taken here. In order to be eonsidered a sueeessful extension of eur- 
rent theory, this iterative approaeh should allow us to prediet observed phenomena. Miller’s 
rule, an empirieal method of attaining higher order eleetrie suseeptibilities in regular non- 
magnetie struetures (i.e., erystals) from the lower eleetrie suseeptibilities, has been shown 
to be effeetive in predieting the nonlinear optieal behavior of these struetures. If valid, 
the approaeh taken here to nonlinear optieal responses should reduee to Miller’s rule in 
the appropriate material regime. As a further test of the predietive power of our approaeh, 
we should be able to show the appearanee of well-known nonlinear optieal effeets sueh as 
frequeney doubling and the Kerr effeet in appropriate media. 
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APPENDIX A: 
Green Functions 


Green functions relevant to the iterative approach of chapter 5 of this work will be derived 
here. 


A.l The Green function for the two dimensional Helmholtz 
equation 

The homogenous Helmholtz equation is given as 


^^f + klf = 0 (A.l) 

where / is the function of interest and ko is a known parameter. The Green function for this 
equation is therefore given as the solution to 

V^G (r, r') +klG (r, /) = 5 {r- r') (A.l) 

Performing a Fourier transform and applying the sifting property of the Dirac delta func¬ 
tion, 



{V^G{r,r')+klG{r,r')) 





(A.3) 



I k^G(r,r') 




ihr' 


From Green’s second identity, 


(A.4) 


J^{fV^g-gV^f)dA 



g(vf-n))di 


(A.5) 
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where D is the domain of the funetions / and g, dD is the bounding set to that domain, and 
n is the normal to the bounding set. Thus, 


/ /vVa= / 

Id Jd JdD \ an an 


(A.6) 


d f d 

where — = V/- H and similar for —. Taking / = and g = G, and furthermore taking 
on on 

D as the entire x-y plane, and assuming that G falls off at least as fast as -, assuming a 

r 

value of zero in the large r limit. 



e^'^'VG(r,r) JA 



(A.7) 



e-^^'^W^G{r,r)dA 



(A.8) 






(A.9) 


The two-dimensional Fourier transform of the Laplaeian of the Green funetion is thus sim¬ 
ply the Fourier transform of the Green funetion multiplied by —k^, where k is the Fourier 
transform parameter. Inserting this result into (A.7), with G^^ I-oc G (r, r') e ^ ''d^r. 


{kl - k^) G (k, r ) = (A. 10) 

The Green funetion may now be attained by inverse Fourier transforming the above expres¬ 
sion. 
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(A.12) 


G 



1 poo „—ik-7' 

^ I 1 _ Jk-r.27 

ilKYJ-ookl-k'^ 


G 



1 |.oc ^ik-ir-r') 
{iTtf J-oo kl-k^ ^ ^ 


(A.13) 


In Cartesian coordinates in k-space (this integral may also be performed in cylindrical 
coordinates), 


G 



-1 

(2;r)2 



^i(^kxPx+kypy^ 

^2 _ p p 

Kq Ky 


dkx dky 


where p = r — r'. 


G 



1 

(2;r)2 




yikxPx 


kx 


{kl-kf) 


dkx dky 


(A. 14) 


(A.15) 


G 




(A.16) 


The denominator of the kx integral has simple poles given by = ± V A:2 _ ^2 jf 
greater than ky, these poles lie on the real kx-axis. If ko is less than ky, the poles lie on the 
imaginary kx-axis. If ko = ky, there is a single pole of second order at the origin. We will 
evaluate (A.16) case by case. 


A.1.1 Case 1: |^o| > \ky\ 

The integrand of the kx integral has poles on the real kx line at kx — — ky. We will 

evaluate it via a contour in the complex kx plane comprised of a large semicircle in the 
upper half of the plane, taking Jk^ — k^ —)■ Jk^—kj + ie. 
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Taking f{z) = 


,KPx 


z+ \lkl-k^ + ie ] [Z- 


k-Q—ky— ie 


, we note that the poles of f{z) 


now lie at ^jk^ — k'^ + ie and —^jk^ — k^ — ie, and that therefore only the simple pole 




kg — ky + ie lies within our contour. From the residue theorem, 


V" V' V 0 > J J J-R JCr 


Ini Res (f(z= J kl - k^ + ie]] = [ f %) dk^ + l_ f (z) dz 

-R 


(A.17) 


where Cr denotes the semicircle in the upper half plane of radius R. 


= / fikx)dk,+ 

J-R 


JzPx 


{z+ xlkl-kj + ie\ (z- Jkl-kj-ie 


dz 


(A.18) 


Taking z = on Cr, 


(■R /-Tt 

f{kx) dkx + 


JRe'^Px 


10 I j^^io + ^ /^2 _ ^2 + ie \ I - Jkl -kj- ie 


iRe^^dd (A. 19) 


ij^^id-\-ipxR{co&{6)-\-i&m{Q)) 


fR rn 

f (kx) dkx + / 

I Re‘^ + xlkl-kj + ie ) iRe ‘0 - Jk^ -kj- ie 


de (A.20) 


= / f{kx)dkx + 
J-R 


iJ^eiQ+iPxRcos,(e) ^-PxRsm{e) 


Re'® + 


Kq 


■ ky A ie 


Re 


iO 


Kq 


kj- 


le 


de (A.21) 


As R —)■ oo, the second integral approaches a value of zero by Jordan’s lemma, since 
sin(0) > 0 on 0 = [0, ;r]. 
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Inserting this result into (A. 22), 
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A.1.2 Case 2: |^o| < |^ 3 ^| 

In this case, f {£) has simple poles at kx — 

large semicircle in the upper half plane, only i 
residue theorem, 

IniRes {j (^z = isjkj- j = j (kx) dkx + f{z)dz (A.30) 

where fiz) = 7 - , , and Cr is again the semicircle of radius R in the 

upper half plane, centered at the origin. Taking z = Re'® on Cr, 


ky—k^. Again closing our contour with a 
Iky—k^ lies within the contour. From the 


■R rlt JRe'^Px 

f{kx)dkx+ 7 —^- , X . —^- iRe'^de (A.31) 

‘ ^ ky — k^ ^Re'® — i^Jky—k^ 


■R pK •D„ie+i/?(cos(e)+!sin(e))pj: 

f{kx)dkx+ — - . -' - - ,-, dO (A.32) 

(^Re'® + i^kj-kl^ (^Re'® - i^k'^-k^ 
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Inserting this result into (A.34), 


(A.37) 


(A.38) 


(A.39) 
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^iPxkx 





(A.40) 


(A.41) 


A.1.3 Case 3: |^o| = |^ 3 ^| 

In this case, our kx integral reduces to 


yikxPx 

kx 


dkx 


(A.42) 


Moving the (now second-order) singularity off of the real axis as for case 1, we consider 


JzPx 


Ic (z + ie) (z-ie) 


rrdz 


(A.43) 


where C is once more the large semicircle in the upper half plane. By the residue theorem, 


Ini Res (/ (z = ze)) = 


^ikxPx 


-—r dkx + 


oizPx 


R {kx + ie) {kx - ie) Jcr {z + ie) {z - ie) 


— dz (A.44) 


where / (z) = 


,lZPx 


Taking z = Re'^ on Cr, 


{z + ie) {z-ie) 


yikxPx 


-—dkx + 




R {kx -h ie) {kx - ie) Jo [Refe) [Re^^ - ie) 


iRe^^ d 9 


(A.45) 


J-R 


JkxPx rn ^ipxR{cos{e)+ism{e)) 

dkx+ . , . , . — r-^de (A.46) 


{kx -h ie) {kx - ie) Jo {Re^^ + ie) [Re^^ - ie) 
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yikxPx 


(kx + ie) (kx-ie) 


rTt ij^^i6+ip.tRcos{e)-p.cRsm{e} 

dkx+ / - ^dO 

Jo (Re^^ + is) [Re‘^ — ie) 


(A.47) 


As previously, the integral over Cr approaches zero as 7? —)■ oo by Jordan’s lemma. 



pikxPx 

_ Ah 

{kx + ie) {kx-ie) " 


Ini Res {f{z 


ie)) 


The pole at z = ie is simple, and therefore 


(A.48) 



pikxPx 

_ Ah 

{kx + ie) {kx-ie) " 


Ini lirn 

Z — 


^iPxZ 

(z ~ 7--TT 

{z + ie){z-ie) 


(A.49) 



pikxPx 

_ 

{kx + ie) {kx-ie) " 


iKi lim 


z — 


APxZ 


[{z + ie) 


(A.50) 



pikxPx 

_ Ah 

{kx + ie) {kx-ie) " 


Ini 


giPxie 

lie 


(A.51) 



pikxPx 

_ Ak 

{kx + ie) {kx-ie) " 


ne 


e 


(A.52) 


This integral therefore does not converge as £ —)■ 0. We conclude that ky = kp is an unphys¬ 
ical value for this system. 


A.1.4 Returning to the Green function 

From (A. 16), taking all allowed values of ky from the previous three subsections. 



(A.53) 
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G(r,r) 


1 f nie’-^yPy1 f 7te‘^yPye P'^'J^y 
- 9 / - . :- dky H-^ / - , :- dky 

(iTt) J\ky\<\ko\ ^kl-kj (2^) 2|fc,|>|*o| 

(A.54) 


1 

( 2 ^)^ 



■ ikyPy+ipx^ kl-kj 



dky 


+ 




(A.55) 


( 2 ^)^ 


l^ol ;i:e *0 {kypy) + i sin {kypy)) 


i -2 _ ^2 


d k\ 


+ 


Ifcol Tti (cos (kypy + pjcJkl-k^] + i sin (kypy + py Jkl -k^ 


2-|fcol 


kl-kj 


dky 


+ 


Tte (cos (kypy) + i sin {kypy)) 


'l^ol 


kj-kl 


dky I (A.56) 


We note now that the first and third integrals above may be combined to yield a single 
integral which is symmetric about the origin, while the second integral is already symmetric 
about the origin. Further splitting these integrals into even and odd components, 
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(A.57) 


where Fi is the entire real ky line exeept for the interval [— |A:o|, |A:o|], whieh is symmetrie 
about the origin. Sinee the seeond and fourth integrals above are odd in ky and are integrated 
symmetrieally with respeet to the origin, they integrate to zero. We are left with 


G(r,/) = 


Ke 


Px\J^ ^ 


eos (kypy) 




kj-kl 


dky “h 


Mo\ Tti eos (kypy + pxsjkl-kf^ 
'-Ifcol 


dkv 


kl-kj 


(A.58) 


The remaining integrals are eaeh even and integrated symmetrieally with respeet to the 
origin. 


G 



2 

(2;r)2 



Ke 


Px\Jk'^ k 


eos {kypyj 


kj-kl 


dky + 


TT/eos [kypy + p. 


kl-kj. 


kl-kj 



(A.59) 

Without loss of generality, let p lie along the ky axis. Then px = 0 and py = p = |r —r'|. 
We are left with 
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/ 


G(r,r) = 


{IkY 


KCOS (kyp) 


ko 


k) -1 


: dky “ 1 “ 


;r/cos (kyp) 


\ 




'■ dky 


— 1^1 

^0 / 


(A.60) 


rC 

Let M = 7 ^. Note that for ky —)■ |A:o|, u —)■ ±1. 
ko 


1 ( r 7tcos{kopu) 7ticos{kopu 

G(r,rj = ^iyi ^ ^ kodu+ / ^ ^ kodu 


ko\/ —I 


lo koVl — 


(A.61) 


Since the second integral above is even in u. 


1 ( r7tcos{kopu) 7ticos{kopu 

G(r,rj = ^ 1^1 ^— ^=^kodu+ I —— , kodu 


ko\/u^ — 1 


'0 koV 1—u^ 


(A.62) 


For the second integral above, let u = cos {t). 


G[Yr) = ^ 


27r I J\ y/iP- — \ 


cos(kopu) , —icos(kopcos(t)) . , i 

^ ^ ^ du+ I - ) ^ ^ sin (t) dt 1 (A.63) 


\l sin^ it) 


^.y 1 f f°°cos(kopu) , , [2 / NX , 

G{r,r) = — { / — ^ du + i cos {kop cos (t)) dt 

2% \J\ a/m2 - 1 Jo 


(A.64) 


The first integral above is an integral representation of the Bessel function of the second 
kind, order 0, 


rcos{kopu) n ^ /A 

/ /9 , du = -No{kop) (A.65) 

Ji Vm^ -1 2 

While the second integral of (A.63) is an integral representation of the Bessel function of 
the first kind, order 0, 
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(A.66) 


/* 2 TT 

/ cos {kop cos (t)) dt = —jQ^kp) 
Jo 2 


Inserting (A.64) and (A.65) into (A.63), 


G (r,/) = ^ (kop) + y7o (kop)^ 
G {r, r) = 2 (-^0 (kop) iNq (kop)) 
G{r,r) = ^-H^^\kop) 


where H, 


( 1 ) 


0 


is the first Hankel funetion of order zero. 


(A.67) 

(A.68) 

(A.69) 
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APPENDIX B: 
MATLAB Simulations 


This appendix will present the MATLAB seripts used to simulate the various TO-derived 
struetures studied in Chapter 5 of this work. 


B.l Cylindrical Cloak in Cartesian Mesh 

%% An Iterative Approach: Cylindrical Structure with Cartesian Mesh 

9 - 9 -- 
0 0 

% Input physical constants 
%clearvars; 

run( 'PhysicalConstants.m' ); 

9 - 9 -- 
o o 

9 - 9 -- 
o o 

% Input known material values 
a = 1; 
b = 1.5; 

R = (b-a)/b; 

9 - 9 -- 
o 0 

% Input meshgrid for the region of interest, 
increment = 0.03; 

X = cat(2,-b:increment:-increment,increment:increment:b); 
y = cat{2,-b:increment:-increment,increment:increment:b); 

[X,Y] = meshgrid(x,y); 

9 - 9 -- 
o o 

% Convert to cylindrical coordinates 

r = sqrt{X.^2+Y.^2) ; 

origin = size(x,2)/2; 

theta = zeros(size(x,2),size(x,2)); 

for n = 1:1:size (x,2)/2 

for m = 1:1:size(x,2)/2 

theta(n,m) = atan(y{n)/x(m))+pi; 

theta(n+size(x,2)/2,m) = atan(y{n+size(x,2)/2)/x(m))+pi; 
theta(n,m+size(x, 2)/2) = atan(y(n)/x(m+size(x,2)/2))+2*pi; 
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theta(n+size(x,2)/2,m+size(x,2)/2) = ... 

atan (y{n+size(x,2)/2)/x(m+size(x, 2)/2) ) ; 


end 


end 


% Input the coordinate transformation from $r$ to $r"' 
r_prime = (r-a) ./R; 
for n = 1 :1: size (x,2) 

for m = 1:1:size(x,2) 
if r(n,m) <= a 

r_prime(n,m) = 0; 

end 

end 

end 


{ ’ }$. 


% Input the transformation optics derived susceptibilities. Assume vacuum 
% outside of the "shield." 

L_ll = zeros (size(x,2),size(x,2)); 

L_12 = zeros(size(x,2),size(x, 2)) ; 

L_13 = zeros (size(x,2),size(x,2)); 

L_22 = zeros(size(x,2),size(x,2)); 

L_23 = zeros (size(x,2),size(x,2)); 

L_33 = ones(size(x,2),size(x,2)); 

L = cell(size(x,2),size(x,2)); 

Chi_ll = zeros(size(x,2),size(x,2)); 

Chi_12 = zeros(size(x,2),size(x,2)); 

Chi_13 = zeros(size(x,2) , size(x, 2)) ; 

Chi_22 = zeros(size(x,2),size(x,2)); 

Chi_21 = zeros(size(x,2),size(x,2)); 

Chi_23 = zeros(size(x,2),size(x,2)); 

Chi_31 = zeros(size(x,2),size(x, 2)) ; 

Chi_32 = zeros (size(x,2),size(x,2)); 

Chi_33 = zeros(size(x,2),size(x,2)); 
for n = 1:1:size (x,2); 

for m = 1:1:size (x,2) 

L_11 (n, m) = R+ (a/ r_prime (n,m) )* (sin (theta (n,m) ) ) '^2; 

L_12(n,m) = -a*cos(theta(n,m))*sin(theta(n,m))/r_prime(n,m); 
L_22 (n, m) = R+ (a/ r_prime (n,m) ) * (cos (theta (n,m) ) ) '^2; 

end 
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end 

for n = 1:1:size (x,2) 

for m = 1:1:size {x,2) 

L{n,m} = [L_11(n,m) , L_12{n,m) ,L_13{n,m);L_12{n,m),L_2 2{n,m), . . . 
L_2 3(n, m);L_13(n,m),L_2 3(n,m),L_33(n,m)]; 

end 


end 


Chi = cell {size (x,2),size(x,2)); 
for n = 1:1:size(x,2) 

for m = 1:1:size(x,2) 
if r(n,m) > b 

Chi{n,m} = zeros(3,3); 
else if r(n,m) <= a 

Chi{n,in} = zeros(3,3); 

else 

Chi{n,in} = (L{n,m} * ( (L{n,m}) 
r (n,in) ) ) -eye (3) ; 

end 

end 

end 


end 

for n = 
for 


1:1:size (x,2) 
m = 1:1: size (. 
Chi_l1(n,m) = 
Chi_12(n,m) = 
Chi_13(n,m) = 
Chi_21(n,m) = 
Chi_22(n,m) = 
Chi_23(n,m) = 
Chi_31(n,m) = 
Chi_32(n,m) = 
Chi_33(n,m) = 


, 2 ) 

Chi{n,m}(1,1) 
Chi{n,m}(1,2) 
Chi{n,m}(1,3) 
Chi{n,m}(2,1) 
Chi{n,m}(2,2) 
Chi{n,m}(2,3) 
Chi{n,m}(3,1) 
Chi{n,m}(3,2) 
Chi{n,m}(3,3) 


end 


))•*(r_prime(n,m)/(R*. .. 


end 

9 - 9 -- 
0 o 

% input the incident fields, as well as the functions $F_1(E_0,H_0)$ and 
% $F_2(H_0,E_0)$. Begin by choosing a value for $\omega_0$. Take $10^{11} Hz$. 
k_0 = le-2; 


9 - 9 -- 
o o 
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123 
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137 

138 
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141 
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144 


% Input the initial (frequency domain) electric field 
E0_z = 2*pi . *exp (li*k_0.9cX) ; 

9 - 9 -- 
o o 

% Compute the tensor contractions appearing in $F_1S and $F_2$ 


9 - 9 -- 
o o 

% Perform the integration indicated in equation (3.40) to obtain the 
% first-order response electric field. Begin on the x-axis outside the 
% shield. That is, let $\vec{r} = -2a\hat{x}$. 
rminusrprime = cell(size(x,2),size (x,2)); 

G1 = cell(size(x,2),size(x,2)); 

51 = 3*X.^2.*Y.^4 + 3*X.^4.*Y.^2-r.^5+X.^6+Y.''6; 

52 = X.^2-2*r.^3+Y.^2; 

Chi_22_x = X.*(6*X.^2.*r-r.^3+2*Y.^2)./Sl-X.*(12*X.^2.*Y.^2-5*r.''3+6*. .. 

X.''4 + 6*Y."4) .* (X. "2+Y. "2-r. ■'5+2 . *X.-^2 . *r . "3+Y.-'4) ./ (Sl.''2) ; 

Chi_12_y = X.*S2./S1-2*X.*Y.-^2.*(3*r-l) ./S1-X.*Y.^2.*S2.*(12*X.^2.*Y.^2-... 
5*r.''3 + 6*X."4 + 6*Y."4) ./(SI.''2); 


for n = 1:1:size(x,2) 

for m = 1:1:size(x,2) 
if r(n,m) > b 

Chi_22_x(n,m) = 0; 

Chi_12_y(n,m) = 0; 
elseif r(n,m) <= a 

Chi_22_x(n,m) = 0; 

Chi_12_y(n,m) = 0; 

end 

end 

end 

for n = 1:1:size (x,2); 

for m = 1:1:size(x,2); 

rminusrprime {n, m} = sqrt ( (X (n, m) -X) .^2+(Y(n,m)-Y) .'^2); 
Gl{n,m} = besselh(0,l,k_0.*rminusrprime{n,m}); 

G1 {n, m} (n, m) = 0 ; 
end 

end 

9 - 9 - 
o o 

% Compute the integral (4.58) 

F = cell(size(x,2),size(x,2)); 
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for n = 1:1: size (x,2) 

for m = 1:1:size (x,2) 

F{n,in} = pi*li/2*Gl{n,m}.*exp( li9rk_0*X) . * {k_0 '^2*Chi_33-li*k_0* . . . 
Chi_2 2_x+k_0^2*Chi_22 + li*k_09rChi_12_y) ; 

end 

end 

E0_z = 2*pi-^exp (li + k_0-^X) ; 

E0_x = 0; 

El_z = zeros (size{x,2),size{x,2)); 
for n = 1:1:size(x,2) 

for m = 1:1:size (x,2) 

El_z(n,m) = increment^2 *trapz(trapz(F{n,m},2)); 

end 

end 

9 - 9 -- 
o o 

% Linear approach. 

E_0 = 1; 

D_x = {Chi_ll+1)*E_0; 

D_y = {Chi_21)*E_0; 

D_r = D_x.*cos(theta)+D_y.*sin(theta); 
for n = 1:1:size(x,2) 

for m = 1:1:size(x,2) 
if r(n,m) <=a 

D_r(n,m) = 0; 

end 

end 

end 

9 - 9 -- 
o o 

% streamline plot the linear field 
streamline(X,Y,D_x,D_y,-b*ones(l,size(x,2)),y); 

9 - 9 -- 
0 0 

% Graphs 

hi = surf(X,Y,Chi_ll); 

title ('$ \chi_{ XX }$ for Cylindrical Redirection Structureinterpreter . 
'latex' , 'fontsize' ,14) ; 

xlabel ( 'X (meters )' , ' interpreter ' , 'latex ' , ' fontsize ' ,12 ); 
ylabel( 'y (meters)' , 'interpreter' , 'latex' , 'fontsize' ,12); 
zlabel( '$\chi_{xx}$' , 'interpreter' , 'latex' , 'fontsize' , 14) ; 
figure 
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surf (X, Y, abs (E0_z+El_z) / (29rpi) ) 

xlabel{ 'X (meters)' , 'interpreter' , 'latex ') ; 

ylabel {'y (meters)','interpreter' , 'latex ') ; 

zlabel( '$\frac{\left|E_0+E_1\right|}{\left|E_0\right|}$' , 'interpreter' , ... 

'latex' , 'fonts!ze' ,14); 
figure 

surf(X,Y,abs(El_z)/(2*pi)) ; 

title ('1st Order Response Field Amplitude Relative to $\left|E_0\right|$' , ... 

'interpreter' , 'latex' ); 
xlabel( 'X (meters)' , 'interpreter' , 'latex ') ; 
ylabel( 'y (meters)' , 'interpreter','latex ') ; 

zlabel( '$\frac{\left|E_1\right| }{\left|E_0\right| }$' , 'interpreter' , . . . 

'latex' , 'fontsize' ,14); 


B.2 Cylindrical Structure in Cylindrical Mesh 

%% An Iterative Approach: Cylindrical Structure in Cylindrical Mesh 

9 - 9 -- 
o o 

% Input physical constants 
clearvars; 

%run('Physicaiconstants.m'); 

9 - 9 -- 
0 0 

% Input known material parameters 
a = 1; 
b = 1.5; 

R = (b-a)/b; 

9 - 9 - 
o o 

% Input cylindrical meshgrid for the region of interest 
rstep = 0.02; 
rmax = 2; 
phistep = pi/50; 

[r,phi] = meshgrid{0:rstep:rmax,0:phistep:2*pi) ; 

9 - 9 - 
o 0 

% Input relative permittivity in cylindrical coordinates. 

eps_rr = ((r-a)./r).^2; 

eps_pp = ones(size(r)); 

eps_zz = (1/R)^2.*ones(size(r)); 
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n = 1: 

1:size (phi,1) 


for m 

= 1:1:size(r,2) 


if 

r(n, m)>b 



eps_rr(n,m) = 

1; 


eps_pp(n,m) = 

1; 


eps_zz(n,m) = 

1; 

el 

seif r (n,m)<=a 



eps_rr(n,m) = 

1; 


eps_pp(n,m) = 

1; 


eps_zz(n,m) = 

1; 


end 

end 


end 

9 - 9 -- 
o o 


% Convert 
eps_xx = 
eps_yy = 
eps_xy = 
chi_xx = 
chi_yy = 
chi_xy = 
chi_yx = 
chi_zz = 

9 - 9 -- 
0 0 


to Cartesian coordinates. Compute the 
eps_rr .* (cos (phi) ) . ''2+eps_pp . * (sin (phi) 
eps_rr .*(sin(phi) ) . ''2+eps_pp . * (cos (phi) 
(eps_rr-eps_pp).*sin(phi).*cos(phi); 
eps_xx-l; 
eps_yy-l; 
eps_xy; 
eps_xy; 
eps_z z-1; 


susceptibility. 
) .'' 2 ; 

) .'' 2 ; 


% Plot the XX component of the permittivity for comparison with the purely 
% Cartesian approach 

hi = surf(r.*cos(phi),r.*sin(phi),chi_xx); 
hold on 

title ('$\chi_{xx}$ for Cylindrical Redirection Structureinterpreter . 
'latex' , 'fontsize' ,14); 

xlabel( 'X (meters)' , 'interpreter' , 'latex' , 'fontsize' ,12); 
ylabel ('y (meters)','interpreter' , 'latex' , 'fontsize' , 12) ; 
zlabel( '$\chi_{xx}$' , 'interpreter' , 'latex' , 'fontsize' ,14); 
hold off 


o g_ 
o o 


% Input incident wavenumber and electric field 
k_0 = le-1; 

E0_z = 2*pi*exp(li*k_0*r.*cos(phi)); 

9 - 9 - 
o o 
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% Input the terms to be integrated 
X = r . 9rcos (phi) ; 

Y = r.*sin(phi); 

51 = 3*X.'^2.*Y.^4 + 3*X.^4.*Y.'^2-r.^5+X.^6+Y.^6; 

52 = X.^2-2*r.^3+Y.^2; 

chi_yy_x = X.* ( 6*X . ^2 . *r-r . ^ 3 + 2 * Y . ^2 ) ./Sl-X.*(12*X.^2.*Y.^2-5*r.'^3 + 6*. . . 

X. "'4 + 6*Y."4) (X. "2+Y. "2-r.'^5 + 2 . *X. "'2 . *r . "3+Y. "'4) ./ (Sl."'2) ; 

chi_xy_y = X . *S2 . / S1-2 *X . * Y . ^^2 . * (3*r-l) ./S1-X.*Y.^2.*S2.*(12*X.^2.*. . . 

Y. ^2-5*r."3+6*X."4+6*Y."4)./(Sl."2); 
for n = 1:1:size (r,1) 

for m = 1:1:size (r,2) 
if r(n,m) > b 

chi_yy_x(n,m) = 0; 
chi_xy_y(n, m) = 0; 
eiseif r(n,m) <= a 

chi_yy_x(n, m) = 0; 
chi_xy_y(n, m) = 0; 

end 

end 

end 

rminusrprime = cell(size(r)); 

G = cell(size(r)); 
rdphi = phistep9rr (1, : ) ; 
for n = 1:1:size (r,1) 

for m = 1:1:size (r,2) 

rminusrprime{n,m} = sqrt((X(n,m)-X).^2+(Y(n,m)-Y).^2); 

G{n,m} = -besselh(0,1,k_0.*rminusrprime{n,m}); 

G{n,m}(n,m) = 0; 

end 

end 

F = cell(size(r)); 
for n = 1:1:size (r,1) 

for m = 1:1:size (r,2) 

F{n,m} = pi*li/29cG{n,m}.*exp (li9rk_0*X) . * (k_0 ^2*chi_zz-li*k_09r . . 
chi_yy_x+k_0 ^2 *chi_yy+1 i*k_0 *chi_xy_y) ; 

end 

end 

El_z = zeros(size(r)); 
for n = 1:1:size (r,1) 
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for m = 1:1: size (r,2) 

El_z(n,m) = rstep*trapz(rdphi.rapz(F{n,m},l)); 

end 

end 

9 - 9 -- 
o o 

% Graphs 
figure; 

h2 = surf(X,Y,abs{El_z+E0_z)/{2*pi)); 
hold on 

title (' 1st Order Response Field Amplitude Relative to $\left|E_0\right|$' , ... 

'interpreter' , 'latex' ); 
xlabel{ 'X (meters)' , 'interpreter' , 'latex ') ; 
ylabel {'y (meters)',' interpreter ',' latex ') ; 

zlabel( '$\frac{\left | E_0+E_1\right| }{\left | E_0\right I }$' , 'interpreter' , . . . 
'latex' , 'fontsize' ,14); 

title ('1st Order Combined Field Amplitude Relative to $\left|E_0\right|$' , ... 

'interpreter' , 'latex' ); 
hold off 

XX = -rmax:rstep:rmax; 

E0_c = cat(2,fliplr(E0_z(51,:)),E0_z(1,:)); 

El_c = cat(2,fliplr(El_z(51,:)),El_z(l,:)); 

E0_c(size(r, 2)+1) = []; 

El_c(size(r,2)+1) = []; 
figure 

plot(xx, abs(E0_c+El_c)/(2*pi)) 
hold on 

T = '1st Order Combined Field Along X-axis Relative to $\left|E_0\right|$' ; 

title(T, 'Interpreter' , 'latex' ); 

xlabel( 'X (meters)' , 'interpreter' , 'latex ') ; 

ylabel( '$\frac{\left | E_0+E_1\right| }{\left | E_0\right I }$' , 'interpreter' , . . . 

'latex' , 'fontsize' ,14); 
hold off 


B.3 Reduced Material Parameter Cylindrical Structure 
in Cylindrical Mesh 

%% An Iterative Approach: Reduced Parameter Cylindrical Structure 
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28 

29 
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o, g, 

0 0 

% Input physical constants 
clearvars; 

%run ( ^ Physicaiconst ant s . in ’ ) ; 

9 - 9 -- 
o o 

% Input known material parameters 
a = 1; 
b = 1.5; 

R = (b-a)/b; 

o g, 
o 0 

% Input cylindrical meshgrid for the region of interest 
rstep = 0.008; 
rmax = 2; 
phistep = pi/50; 

[r,phi] = meshgrid{0:rstep:rmax,0:phistep:2*pi) ; 

9 - 9 - 
o o 

% Input relative permittivity in cylindrical coordinates. 

eps_rr = ((r-a)./r).''2; 

eps_pp = ones (size (r)); 

eps_zz = (1/R)^2.*ones(size(r)); 


siz 

e(phi. 

1) 


1:1 

:size ( 

r, 2) 


' (n. 

m) >b 



eps. 

_rr(n. 

m) = 

1; 

eps. 

_PP(n. 

m) = 

1; 

eps. 

_zz(n. 

m) = 

1; 

ill 

r (n, m) 

<=a 


eps. 

_rr(n. 

m) = 

1; 

eps. 

_PP(n. 

m) = 

1; 

eps. 

_zz(n. 

m) = 

1; 


end 

end 

end 

g, g, 

0 0 

% Convert to Cartesian coordinates. Compute the susceptibility. 
eps_xx = eps_rr . * (cos (phi) ) .''2+eps_pp . * (sin (phi) ) .''2 ; 
eps_yy = eps_rr . * (sin (phi) ) .''2+eps_pp . * (cos (phi) ) .''2 ; 
eps_xy = (eps_rr-eps_pp).*sin(phi).*cos(phi); 
chi_xx = eps_xx-l; 
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chi_yy = eps_yy-l; 
chi_xy = eps_xy; 
chi_yx = eps_xy; 
chi_zz = eps_zz-l; 

9 - 9 -- 
o o 

% Plot the XX component of the permittivity for comparison with the purely 
% Cartesian approach 

hi = surf(r.*cos(phi),r.*sin(phi),chi_xx); 
hold on 

title ('$\chi_{XX}$ for Cylindrical Redirection Structureinterpreter . 
'latex' , 'fontsize' ,14); 

xlabel{ 'X (meters) ' , 'interpreter' , 'latex' , 'fontsize' ,12); 
ylabel ('y (meters)','interpreter' , 'latex' , 'fontsize' ,12); 
zlabel( '$\chi_{xx}$' , 'interpreter' , 'latex' , 'fontsize' , 14) ; 
hold off 

9 - 9 -- 
o o 

% Input incident wavenumber and electric field 
k_0 = le-1; 

E0_z = 2*pi*exp(li*k_0*r.*cos(phi)); 

9 - 9 -- 
0 o 

% Input the terms to be integrated 
X = r.*cos(phi); 

Y = r.*sin(phi); 

chi_yy_x = (2*X.*Y.^2) .* (3*r-2) ./ (r.^6); 

chi_xy_y = -X . * (2 *X .''2 . *r-4 *Y . ^2 . *r-X . ^2 + 3*Y . ^2 ) . / (r . 6) ; 
for n = 1:1:size(r,1) 

for m = 1:1:size (r,2) 
if r(n,m) > b 

chi_yy_x(n,m) = 0; 
chi_xy_y(n,m) = 0; 
elseif r(n,m) <= a 

chi_yy_x(n,m) = 0; 
chi_xy_y(n,m) = 0; 

end 

end 

end 

rminusrprime = cell(size(r)); 

G = cell(size(r)); 
rdphi = phistep*r(1, :) ; 
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for n = 1:1:size (r,1) 

for m = 1:1: size (r,2) 

rminusrprime{n,m} = sqrt((X(n,m)-X).^2+(Y{n,m)-Y).^2); 

G{n,m} = -besselh(0,1,k_0.*rminusrprime{n,m}); 

G{n,m}(n,m) = 0; 

end 

end 

F = cell{size(r)); 
for n = 1:1:size (r,1) 

for m = 1:1:size (r,2) 

F{n,m} = pi*li/29cG{n,m}.*exp (li9rk_0*X) . * (k_0 ^2*chi_zz-li*k_09r . . . 
chi_yy_x+k_0 ^ 2 chi_yy +1 i *k_0 *chi_xy_y) ; 

end 

end 

El_z = zeros(size{r)); 
for n = 1:1:size (r,1) 

for m = 1:1: size (r,2) 

El_z(n,m) = rstep^trapz{rdphi.rapz(F{n,m},l)); 

end 

end 

9 - 9 - 
0 0 

% Graphs 
figure; 

h2 = surf(X,Y,abs{El_z+E0_z)/(2*pi)) ; 
hold on 

title ('1st Order Combined Field Amplitude Relative to $\left|E_0\right|$' , ... 

' interpreter' , 'latex' ); 
xlabel{ 'X (meters)' , 'interpreter' , 'latex ') ; 
ylabel( 'y (meters)' , 'interpreter' , 'latex ') ; 

zlabel( '$\frac{\left|E_0+E_1\right|}{\left|E_0\right|}$' , 'interpreter' , ... 
'latex' , 'fontsize' ,14); 

title (' 1st Order Combined Field Amplitude Relative to $\left|E_0\right|$' , ... 

'interpreter' , 'latex' ); 
hold off 

XX = -rmax:rstep:rmax; 

E0_c = cat(2,fliplr(E0_z(51,:)),E0_z(1,:)); 

El_c = cat(2,fliplr(El_z(51,:)),El_z(l,:)); 

E0_c(size(r, 2)+1) = []; 

El_c(size(r, 2)+1) = []; 
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figure 

plot(xx,abs(E0_c+El_c)/(2*pi)) 
hold on 

T = ’1st Order Combined Field Along X-axis Relative to $\left|E_0\right|$' ; 

title (T, 'Interpreter' , 'latex ') ; 

xlabel{ 'X {meters)' , 'interpreter' , 'latex ') ; 

ylabel( '$\frac{\left|E_0+E_1\right| }{\left|E_0\right I }$ ' , ’interpreter' , . . . 

'latex' , 'fonts!ze' ,14); 
hold off 


B.4 Square Cloak 


%% An Iterative Approach: Symbolic Method for Square Cloak 

o, g, 

0 o 

% Input physical constants 
clearvars; 

run( 'PhysicaiconStants.m' ); 

9 - 9 -- 
o o 

% Input predetermined structure shape. Let $S_1$ and SS_2$ be the inner and 
% outer half side lengths, respectively, 
si = 1; 
s 2 = 2; 

9 - 9 - 
0 o 

% Input symbolic expression for the electric permittivity, as per Pendry et 
% al . Rotate the symbolic expression through $\frac{\pi}{2}$ , $\pi$, and 
% $\frac{3\pi}{2}$ in order to obtain the permittivity throughout the 
% structure, 
syms X y z; 
a = s2/{s2-sl); 

b = symfun (a*y*sl/ {x''2) , [x y] ) ; 
c = symfun(a*(1-sl/x),[x y]); 

epsl = symfun([c/a,-b/a,0;-b/a,(a^2+b^2)/(a*c),0;0,0,a*c],[x y]); 

9 - 9 - 
0 0 

% Input meshgrid for the region of interest, 
increment = 0.03; 

[X,Y] = meshgrid(-1.5*s2:increment:1.5*s2,-1.5*s2:increment:1.5*s2); 

9 - 9 - 
0 0 
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% Determine numeric values of the permittivity over one quarter of the 
% structure, and rotate, as per [8]. 

A_90 = [0,-1,0; 1, 0, 0; 0,0,1] ; 

A_n90 = [0,1,0;-1,0,0;0,0,1]; 

A_180 = [-1,0,0;0,-1,0;0,0,1]; 

A_nl80 = A_180; 

A_270 = [0,1, 0;-l, 0, 0; 0,0, 1] ; 

A_n270 = [0,-1,0;1,0,0;0,0,1]; 
eps2 = symfun(A_90+ ... 

epsl{dot(A_n90*[x;y;0],[1;0;0]),dot(A_n90*[x;y;0],[0;l;0]))*... 

A_90.’,[X y]); 
eps3 = symfun(A_180* ... 

epsl(dot(A_nl80*[x;y;0],[1;0;0]),dot(A_nl80*[x;y;0],[0;1;0]))*A_180 
[ X y ] ) ; 

eps4 = symfun(A_270* ... 

epsl(dot(A_n27 0*[x;y;0], [1;0;0] ),dot{A_n270*[x;y;0], [0;1;0] ) )*A_2 7 0.',... 
[ X y ] ) ; 

chil = symfun(epsl(x,y)-eye(3),[x y]); 
chi2 = symfun(eps2(x,y)-eye{3) , [x y]) ; 
chi3 = symfun(eps3(x,y)-eye(3) , [x y] ) ; 
chi4 = symfun(eps4(x,y)-eye(3) , [x y]) ; 

CHI = cell(size(X,2),size(X,2)); 

CH1_11 = zeros(size(X,2),size(X,2)); 

CHI_12 = zeros(size(X,2),size{X,2)); 

CHI_22 = zeros(size(X,2),size{X,2)); 

CHI_33 = zeros(size(X,2) , size(X, 2)) ; 

G = cell{size(X,2),size(X,2)); 

k_0 = 2e-l; 

for n = 1:1: size (X,2) 

for m = 1:1:size (X,2) 

if atan2(Y(n,m),X(n,m)) >= -pi/4 && atan2(Y(n,m),X(n,m)) < pi/4 

if abs(X(n,m)) == si 

CHI{n,m} = zeros(3,3); 
elseif abs(X(n,m)) < si && abs(Y(n,m)) < si 
CHI{n,m} = zeros{3,3); 
elseif abs(X(n,m)) > s2 || abs(Y(n,m)) > s2 
CHl{n,m} = zeros{3,3); 

else 

CHI{n,m} = double(chil(X(n,m),Y(n,m))); 
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end 


end 

if atan2(Y{n,m),X(n,m)) >= pi/4 && atan2(Y{n,m),X(n,m)) < 3*pi/4 
if abs{Y{n,in)) == si 

CHI{n,in} = zeros(3,3); 
elseif abs(X(n,m)) < si && abs{Y(n,m)) < si 

CHI{n,in} = zeros(3,3); 
elseif abs(X(n,m)) > s2 || abs{Y(n,m)) > s2 
CHI{n,in} = zeros(3,3); 


else 


CHI {n,in} = 

double (chi2 {X{n,in) ,Y (n,m) ) ) ; 

end 


end 



if atan2(Y{n,m),X(n,m)) >= 3*pi/4 || atan2{Y{n,in),X{n,m)) < -3*pi/4 
if abs{X{n,in)) == si 

CHI{n,in} = zeros{3,3); 
elseif abs(X(n,m)) < si && abs{Y(n,m)) < si 

CHI{n,in} = zeros(3,3); 
elseif abs(X(n,m)) > s2 || abs{Y(n,m)) > s2 
CHI{n,in} = zeros(3,3); 


else 


CHI {n,in} = 

double (chi3 {X{n,in) ,Y (n,m) ) ) ; 

end 


end 



if atan2(Y{n,m),X(n,m)) >= -3*pi/4 && atan2{Y{n,m),X{n,m)) < -pi/4 
if abs{Y{n,in)) == si 

CHI{n,in} = zeros{3,3); 
elseif abs(X(n,m)) < si && abs{Y(n,m)) < si 

CHl{n,in} = zeros{3,3); 
elseif abs(X(n,m)) > s2 || abs{Y(n,m)) > s2 
CHI{n,in} = zeros(3,3); 


else 


CHI {n,in} = 

double (chi4 {X{n,in) ,Y (n,m) ) ) ; 

end 



end 

CHI_ll(n,m) = CHl{n,m} {1,1); 
CHI_12(n,m) = CHI{n,m} {1,2); 
CHI_22(n,m) = CHI{n,m} {2,2); 
CHI_33(n,m) = CHI{n,m} {3,3); 
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G{n,in} = besselh (0, 1, k_0.9csqrt { {X (n, m)-X) . ^2+(Y (n, m)-Y) . ^2 ) ) ; 
G{n,m} (n,m) = 0; 


end 


end 


l;0] 
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% Compute the derivatives of the susceptibility appearing in (4.60) 
chil_12 = symfun(conj(dot(chil(x, y)*[0; 1; 0 ] , [1; i 
chil_12_y = symfun(diff(chil_12(x,y),y),[x y]); 
chi2_12 = symfun(dot(chi2(x,y)*[0;l;0],[1;0;0]), 
chi2_12_y = symfun(diff(chi2_12(x,y),y), [x y]) ; 
chi3_12 = symfun(dot(chi3(x,y)*[0;l;0], [1;0;0]), 
chi3_12_y = symfun(diff(chi3_12(x,y),y),[x y]); 
chi4_12 = symfun(dot(chi4(x,y)*[0;l;0], [1;0;0]), 
chi4_12_y = symfun(diff(chi4_12(x,y),y),[x y]); 
chil_2 2 = symfun(conj(dot(chil(x,y)*[0;1;0] , [Op 
chil_22_x = symfun(diff(chil_22(x,y),x),[x y]); 
chi2_22 = symfun(dot(chi2(x,y)*[0;l;0],[0;1;0]), 
chi2_22_x = symfun(diff(chi2_22(x,y),x), [x y]) ; 
chi3_2 2 = symfun(dot(chi3(x,y)*[0;l;0], [0;1;0]), 
chi3_22_x = symfun(diff(chi3_22(x,y),x),[x y]); 
chi4_22 = symfun(dot(chi4(x,y)*[0;1;0], 
chi4_22_x = symfun(diff(chi4_22(x,y),x) 

CHI_12_y = zeros(size(X,2),size(X,2)); 

CHI_22_x = zeros(size(X,2),size(X,2)); 
for n = 1:1:size (X,2) 

for m = 1:1:size(X,2) 

if atan2(Y(n,m),X(n,m)) >= -pi/4 && atan2(Y(n,m),X(n,m)) < pi/4 

if abs(X(n,m)) == si 
CHI_12_y(n,m) = 0; 

CHI_22_x(n,m) = 0; 

elseif abs(X(n,m)) < si && abs(Y(n,m)) < si 
CHI_12_y(n,m) = 0; 

CHl_22_x(n,m) = 0; 

elseif abs(X(n,m)) > s2 || abs(Y(n,m)) > s2 
CHl_12_y(n,m) = 0; 

CHI_22_x(n,m) = 0; 

else 

CHI_12_y (n, m) = double (chil_12_y (X(n,m) ,Y (n,m) ) ) ; 
CHI_22_x(n,m) = double(chil_22_x(X(n,m),Y(n,m))); 
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end 


end 

if atan2(Y{n,m),X(n,m)) >= pi/4 && atan2(Y{n,m),X(n,m)) < 3*pi/4 
if abs{X{n,in)) == si 
CHI_12_y(n,m) = 0; 

CHI_22_x(n,m) = 0; 

elseif abs(X(n,m)) < si && abs(Y(n,m)) < si 
CHI_12_y(n,m) = 0; 

CHI_22_x(n,m) = 0; 

elseif abs(X(n,m)) > s2 || abs(Y(n,m)) > s2 
CHI_12_y(n,m) = 0; 

CHI_22_x(n,m) = 0; 

else 

CHI_12_y (n, m) = double {chi2_12_y (X{n,m) ,Y (n,m) ) ) ; 

CHI_22_x(n,m) = double{chi2_22_x(X{n,m),Y(n,m))); 

end 

end 

if atan2(Y{n,in),X(n,m)) >=39rpi/4 || atan2{Y(n,in),X(n,m)) < -39rpi/4 

if abs{X{n,in)) == si 
CHI_12_y(n,m) = 0; 

CHI_22_x(n,m) = 0; 

elseif abs(X(n,m)) < si && abs(Y(n,m)) < si 
CHI_12_y(n,m) = 0; 

CHl_22_x(n,m) = 0; 

elseif abs(X(n,m)) > s2 || abs{Y(n,m)) > s2 
CHI_12_y(n,m) = 0; 

CHI_22_x(n,m) = 0; 

else 

CHI_12_y(n,m) = double{chi3_12_y(X{n,m),Y(n,m))); 

CHI_22_x(n,m) = double{chi3_22_x(X{n,m),Y(n,m))); 

end 

end 

if atan2 (Y {n, m) , X (n, m) ) >= -39rpi/4 && atan2 (Y (n, m) , X {n, m) ) < -pi/4 
if abs{X(n,m)) == si 
CHI_12_y(n,m) = 0; 

CHI_22_x(n,m) = 0; 

elseif abs(X(n,m)) < si && abs{Y(n,m)) < si 
CHI_12_y(n,m) = 0; 

CHI_22_x(n,m) = 0; 
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elseif abs(X(n,m) 
CHI_12_y (n,m) 
CHI_22_x (n, m) 

else 

CHI_12_y (n,m) 
CHI_22_x (n, m) 

end 

end 

end 


> s2 I I abs(Y(n,m)) > s2 
= 0 ; 

= 0 ; 

= double {chi4_12_y (X{n,m),Y(n,m) 
= double {chi4_22_x (X{n,m),Y(n,m) 


end 

F = cell{si2e(X,2),size(X,2)); 
for n = 1:1:size (X,2) 

for m = 1:1:size (X,2) 

F{n,in} = pi*li/2*G{n,m}.*exp(li*k_0*X) .*(k_0 ^2*CHl_33-li*k_0* . . . 
CHl_2 2_x+k_0^2*CHI_22 + li*k_0*CHl_12_y) ; 

end 


end 

El_z = zeros(size{X,2),size {X,2)); 
for n = 1:1:size (X,2); 

for m = 1:1:size(X,2) 

El_z(n,m) = increment^2 *trapz(trapz(F{n,m},2)); 

end 


end 

E0_z = 29rpi9cexp (li*k_0*X) ; 

9 - 9 -- 
o o 


% Graphs 
surf(X,Y, 

title (' $\chi_{xx}$ for Square Structureinterpreterlatex ') ; 
xlabel { 'X (meters )' , ' interpreter ' , 'latex ') ; 
ylabel{ 'x (meters)' , 'interpreter' , 'latex ') ; 

zlabel( '$\chi_{xx}$' , 'interpreter' , 'latex ' , 'fontsize' ,14); 
figure 

surf(X, Y,abs(E0_z+El_z)/(2*pi)) 

title (' 1st Order Combined Field Amplitude Relative to $\left|E_0\right|$' , ... 

'interpreter' , 'latex' ); 
xlabel( 'X (meters)' , 'interpreter' , 'latex ') ; 
ylabel ('y (meters)',' interpreter' , 'latex ') ; 

zlabel( '$\frac{\left | E_0+E_1\right|}{\left|E_0\right|}$' , 'interpreter' , ... 
'latex' , 'fontsize' ,14); 
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221 figure 

222 plot (X{101, : ) ,abs {E0_z (101, : ) +El_z (101, : ) ) / (2*pi) ) ; 

223 T = '1st Order Combined Field Along X-axis Relative to $\left | E_0\right | $ ' ; 

224 title (T, ' interpreter ' , ' latex ' ) ; 

225 xlabel('x (meters )',' interpreterlatex ') ; 

226 y label ( ' $\frac{\left | E_0+E_1\ right | } { \ left | E_0\ right | } $ ' , 'interpreter ' , . . . 

227 'latex ', ' font size ' , 14 ) ; 
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